1
我想將matlab代碼轉換爲C.爲此，我將C代碼與Intel MKL庫連接起來，幷包含「mkl_lapacke.h」。 Matlab代碼包含：Lapack和Matlab之間的不同結果
>>A=mldivide(A1,A2)其中A1和A2都是正方形10x10實矩陣。 這可以被解釋爲系統A1 * X = A2
在C代碼，我打電話Dgesv如下的溶液： info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, n, nrhs, a, lda, ipiv,b, ldb); 其中LDA = 10，N = 10和nrhs = 10
問題是由Matlab和Lapack返回的10x10解決方案非常不同！ 這裏是A1 =一個代碼和A2 = B
#include <stdlib.h>
#include <stdio.h>
#include "mkl_lapacke.h"
/* Auxiliary routines prototypes */
extern void print_matrix(char* desc, MKL_INT m, MKL_INT n, double* a, MKL_INT lda);
extern void print_int_vector(char* desc, MKL_INT n, MKL_INT* a);
/* Parameters */
#define N 10
#define NRHS 10
#define LDA N
#define LDB NRHS
/* Main program */
int main() {
/* Locals */
MKL_INT n = N, nrhs = NRHS, lda = LDA, ldb = LDB, info;
/* Local arrays */
MKL_INT ipiv[N];
double a[LDA*N] = {
-0.0091, 0.1024, -0.2640, -0.0956, 0.0635, -0.1776, 0.1257, 0.1048, -0.0869, 0.0106,
-0.0865, 0.2401, 0.0455, -0.0483, -0.2640, 0.3985, 0.1095, -0.2429, 0.1452, -0.0629,
-0.0428, 0.1669, -0.0239, -0.0877, -0.0893, 0.2085, -0.2095, -0.0423, 0.0712, 0.0051,
-0.0458, 0.0043, 0.3219, 0.1583, -0.1277, -0.0598, 0.2033, -0.1075, -0.0131, -0.0277,
-0.0597, 0.2190, 0.0053, 0.0084, -0.0741, -0.0993, 0.3048, -0.0046, -0.0718, -0.0055,
0.0538, -0.0734, -0.2116, -0.0733, 0.0203, 0.2163, 0.0991, -0.1309, 0.1299, -0.0564,
-0.0415, 0.1569, -0.0053, -0.0754, -0.0855, 0.1912, -0.2020, -0.0347, 0.0524, 0.0122,
0.0648, -0.1265, -0.1628, -0.0357, 0.0592, 0.1129, 0.0953, -0.0884, 0.0892, -0.0431,
0.0446, -0.2029, 0.1323, 0.0604, 0.0271, 0.1125, -0.1788, -0.0454, 0.0663, -0.0126,
0.0241, -0.1181, 0.1255, 0.0281, -0.0157, 0.1600, -0.2448, -0.0524, 0.0855, 0.0092,};
double b[LDB*N] = {
-0.2225, -0.2789, 0.1338, -0.3709, -0.4954, -0.1445, -0.0116, 0.0254, 0.0118, 0.0098,
0.0362, -0.3659, -0.1204, -0.0500, 0.1276, -0.0473, -0.2388, 0.0701, -0.3668, -0.0480,
0.2351, 0.0922, -0.0670, -0.1074, 0.2423, -0.3811, 0.0791, -0.2176, -0.0391, 0.0532,
-0.0023, -0.2109, 0.0767, -0.1575, 0.2569, -0.1005, 0.2427, 0.3022, 0.0923, -0.0445,
0.4103, 0.3612, 0.0651, -0.0481, 0.1001, 0.5006, -0.1107, 0.3178, -0.0713, 0.4568,
0.1862, -0.3224, 0.0601, 0.1015, -0.2129, 0.0320, -0.1459, -0.0723, 0.3412, 0.0431,
0.1613, 0.3168, 0.0876, -0.0442, -0.2465, -0.1598, -0.1102, 0.2010, 0.0080, -0.0619,
0.0929, 0.1286, -0.2801, 0.0119, -0.1908, 0.0509, 0.2731, 0.1054, -0.1830, 0.0112,
-0.1971, -0.1049, -0.0354, 0.5010, 0.0685, -0.2606, 0.0225, 0.0164, -0.0140, -0.0002,
0.0452, -0.2061, 0.2058, 0.0156, 0.0198, -0.0294, 0.0453, -0.1110, 0.0098, 0.0145,
};
/* Executable statements */
printf("LAPACKE_dgesv (row-major, high-level) Example Program Results\n");
/* Solve the equations A*X = B */
info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, n, nrhs, a, lda, ipiv,
b, ldb);
/* Check for the exact singularity */
if(info > 0) {
printf("The diagonal element of the triangular factor of A,\n");
printf("U(%i,%i) is zero, so that A is singular;\n", info, info);
printf("the solution could not be computed.\n");
exit(1);
}
/* Print solution */
print_matrix("Solution", n, nrhs, b, ldb);
/* Print details of LU factorization */
print_matrix("Details of LU factorization", n, n, a, lda);
/* Print pivot indices */
print_int_vector("Pivot indices", n, ipiv);
exit(0);
} /* End of LAPACKE_dgesv Example */
/* Auxiliary routine: printing a matrix */
void print_matrix(char* desc, MKL_INT m, MKL_INT n, double* a, MKL_INT lda) {
MKL_INT i, j;
printf("\n %s\n", desc);
for(i = 0; i < m; i++) {
for(j = 0; j < n; j++) printf(" %6.2f", a[i*lda+j]);
printf("\n");
}
}
/* Auxiliary routine: printing a vector of integers */
void print_int_vector(char* desc, MKL_INT n, MKL_INT* a) {
MKL_INT j;
printf("\n %s\n", desc);
for(j = 0; j < n; j++) printf(" %6i", a[j]);
printf("\n");
}
通過dgesv返回的解決辦法是：
LAPACKE_dgesv (row-major)
1.0e+03 *
0.3270 -0.5215 0.0049 0.0619 -0.0199 -0.1558 2.9911 1.1247 -5.4283 5.2655
0.0751 -0.2225 0.1936 0.0490 -0.0678 -0.0201 0.2473 0.1422 -0.4608 0.7307
-0.0683 0.3846 -0.4393 -0.0885 0.1620 0.0024 0.2210 -0.0303 -0.3558 -0.2766
0.1779 -0.9302 1.0602 0.2237 -0.3761 -0.0056 -0.4986 0.0816 0.7990 0.7407
-0.1549 0.3202 -0.0615 -0.0345 0.0257 0.0775 -1.3939 -0.5276 2.5444 -2.5409
-0.0069 -0.0202 0.0594 0.0175 -0.0235 0.0140 -0.1871 -0.0481 0.3191 -0.2247
-0.0360 0.1518 -0.1521 -0.0332 0.0600 0.0122 -0.0471 -0.0597 0.0958 -0.3078
0.0675 -0.1360 0.0075 0.0220 0.0272 -0.0318 0.5894 0.1936 -1.0823 1.0735
-0.0129 0.0052 0.0142 -0.0096 0.0355 0.0096 -0.2460 -0.1281 0.4625 -0.4091
-0.0963 0.1961 -0.0244 -0.0417 -0.0032 0.0743 -0.8836 -0.3268 1.5972 -1.6047
而利用Matlab返回的解決辦法是：
1.0e+03 *
0.1224 -0.0783 -0.1534 -0.0092 0.0609 -0.0555 1.3240 0.4477 -2.3813 2.0963
0.0528 -0.1725 0.1739 0.0410 -0.0549 -0.0089 0.0615 0.0637 -0.1206 0.3758
-0.0706 0.3868 -0.4366 -0.0889 0.1543 0.0029 0.2127 -0.0271 -0.3417 -0.2905
0.1813 -0.9290 1.0499 0.2236 -0.3556 -0.0058 -0.4944 0.0666 0.7945 0.7407
-0.0576 0.1085 0.0151 -0.0005 -0.0141 0.0297 -0.6005 -0.2044 1.0941 -1.0311
0.0037 -0.0425 0.0664 0.0211 -0.0263 0.0088 -0.1006 -0.0141 0.1613 -0.0616
-0.0286 0.1352 -0.1456 -0.0306 0.0543 0.0082 0.0190 -0.0308 -0.0253 -0.1830
0.0264 -0.0438 -0.0292 0.0067 0.0455 -0.0119 0.2628 0.0591 -0.4845 0.4438
0.0037 -0.0281 0.0227 -0.0046 0.0308 0.0012 -0.1030 -0.0717 0.2019 -0.1450
-0.0352 0.0629 0.0242 -0.0202 -0.0285 0.0443 -0.3862 -0.1238 0.6877 -0.6572
答案如何不同？你能編輯你的文章以包含輸入和輸出矩陣嗎？你確定你沒有在一個程序中使用你想要的矩陣的轉置嗎？ – eigenchris 2015-02-24 22:07:45
沒有看到你的矩陣和兩者產生的輸出，它是不可能回答你的問題。 – 2015-02-24 22:18:58
我剛剛添加了代碼和結果 – Didon 2015-02-24 22:38:28