最近要使用 dtrsm 进行解方程,试用了一两个例子,但总是有误,归根结底还是对fortran的一些性质不熟悉,这里做一下记录
1、 XTRSM(SIDE, UPLD, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
主要用来解方程,而且花样很多
op( A )*X = alpha*B, or X*op( A ) = alpha*B,
op( A ) = A or op( A ) = A**T.
我使用了别人x学习BLAS库 -- DTRSM_cocoonyang的专栏-CSDN博客
写的一个例子,错误颇多,这里做一些订正
2、fortran中数组 是按列排列的, 二维在内存也是按列排列的,所以 它经常用的 LDA or LDB 就是说的 行数,或者说一段内存按列 走排列,碰到 LDX个就立刻卡断弯折
所以 数组的大小是 A(LDA,列数),B(LDB,列数)
给数组赋值的时候:
a = (double *) malloc(lda * m * sizeof(double));// a size lda *m
b = (double *) malloc(ldb * n * sizeof(double));// b size ldb *n
for (i = 0; i < lda; i++)//lda =4
{
for (j = 0; j <m ; j++)//m =4
{
a[i + j * lda] = xxx; // 列优先,如果是行优先,这个地方变成 a[i*lda + j]
}
}
3 、具体的里一个例子:
void dtrsm_test()
{
double *a;
double alpha;
double *b;
char diag;
int i;
int j;
int lda;
int ldb;
int m;
int n;
char side;
char transa;
char transb;
char uplo;
char functionName[] = "DTRSM";
printf("================================================ \n");
printf("Testing BLAS library function -- %s \n", functionName);
printf("\n");
//printf ( "TEST DTRSM\n" );
printf(" DTRSM solves a linear system involving a triangular\n");
printf(" matrix A and a rectangular matrix B.\n");
printf("\n");
printf(" 1: Solve A * X = alpha * B;\n");
printf(" 2: Solve A' * X = alpha * B;\n");
printf(" 3: Solve X * A = alpha * B;\n");
printf(" 4: Solve X * A' = alpha * B;\n");
/*
Solve A * X = alpha * B.
*/
side = 'L';
uplo = 'U';
transa = 'N';
diag = 'N';
m = 4;
n = 5;
alpha = 2.0;
lda = m;
ldb = m;
//全部按列排列!
// a(lda,m) = 4 *4 or a(m,m)
// b(ldb,n) =4 *5 or b(m,n)
a = (double *) malloc(lda * m * sizeof(double));
b = (double *) malloc(ldb * n * sizeof(double));
for (i = 0; i < lda; i++)//lda =4
{
for (j = 0; j <= i; j++)//m =4
{
a[i + j * lda] = 0.0;
}
for (j = i ; j < m; j++)
{
a[i + j * lda] = (double) (i + j + 2);
}
}
transb = 'N';
for ( i = 0; i < ldb; i++ )//ldb=4
{
for ( j = 0; j <n; j++ )//n=5
{
b[i+j*ldb] = 10 * (i+1) + (j+1);
}
}
printf("Test matrix \n");
r8mat_print(m, m, a, " A = ");
r8mat_print(m, n, b, " B:");
dtrsm_(&side, &uplo, &transa, &diag, &m, &n, &alpha, a, &lda, b, &ldb);
r8mat_print(m, n, b, " X = inv ( A ) * alpha * B:");
free(a);
free(b);
}
4、DTRSM 中 单元三角阵的研究
XTRSM(SIDE, UPLD, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
其他参数其实说的都很明白,只有这个
DIAG
DIAG is CHARACTER*1
On entry, DIAG specifies whether or not A is unit triangular
as follows:
DIAG = 'U' or 'u' A is assumed to be unit triangular.
DIAG = 'N' or 'n' A is not assumed to be unit
triangular.
试验了一下,如果 DIAG 为U,就是单元三角阵,那么A的主对角线都将设为1
你比如A 原本是
但如果你 设定 DIAG 为 u,那么A会自动变成下面,不管你A的主对角线初值是什么