格式化 LAPACKE 的带状矩阵
Format banded matrix for LAPACKE
我正在尝试使用 LAPACK 的 C 接口求解一般带状矩阵,在 Intel 的 MKL 中称为 LAPACKE。我要调用的函数是 *gbsv
,其中 *
表示格式。不幸的是,我发现 非常 很难找到有关如何使用 C 接口格式化带状矩阵的工作示例。如果有人可以为那里的所有 C 用户提供一个工作示例,我向你保证它会有所帮助。
fortran 布局作为示例给出 here,但我不确定如何将其格式化以输入到 LAPACKE。我还应该注意,在我的问题中,我必须即时构建带状矩阵。所以我有 5 个系数,A、B、C、D、E 用于每个 i 节点,必须将它们放入带状矩阵形式,然后传递给 LAPACKE。
函数LAPACKE_dgbsv()
的原型如下:
lapack_int LAPACKE_dgbsv( int matrix_layout, lapack_int n, lapack_int kl,
lapack_int ku, lapack_int nrhs, double* ab,
lapack_int ldab, lapack_int* ipiv, double* b,
lapack_int ldb )
与Lapack的函数dgbsv()
的主要区别在于参数matrix_layout
,可以是LAPACK_ROW_MAJOR
(C序)或LAPACK_COL_MAJOR
(Fortran序)。如果 LAPACK_ROW_MAJOR
,LAPACKE_dgbsv
将转置矩阵,调用 dgbsv()
然后将矩阵转置回 C 排序。
其他参数的含义与函数dgbsv()
相同。如果使用 LAPACK_ROW_MAJOR
,则 dgbsv()
的正确 ldab
将由 LAPACKE_dgbsv()
计算,参数 ldab
可以设置为 n
。但是,就像 dgbsv()
一样,必须为矩阵 ab
分配额外的 space 来存储因式分解的细节。
下面的例子利用LAPACKE_dgbsv()
通过中心有限差分求解一维稳态扩散。考虑了零温度边界条件,并将正弦波之一用作源项以检查正确性。以下程序由gcc main3.c -o main3 -llapacke -llapack -lblas -Wall
编译:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <lapacke.h>
int main(void){
srand (time(NULL));
//size of the matrix
int n=10;
// number of right-hand size
int nrhs=4;
int ku=2;
int kl=2;
// ldab is larger than the number of bands,
// to store the details of factorization
int ldab = 2*kl+ku+1;
//memory initialization
double *a=malloc(n*ldab*sizeof(double));
if(a==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
double *b=malloc(n*nrhs*sizeof(double));
if(b==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
int *ipiv=malloc(n*sizeof(int));
if(ipiv==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
int i,j;
double fact=1*((n+1.)*(n+1.));
//matrix initialization : the different bands
// are stored in rows kl <= j< 2kl+ku+1
for(i=0;i<n;i++){
a[(0+kl)*n+i]=0;
a[(1+kl)*n+i]=-1*fact;
a[(2+kl)*n+i]=2*fact;
a[(3+kl)*n+i]=-1*fact;
a[(4+kl)*n+i]=0;
//initialize source terms
for(j=0;j<nrhs;j++){
b[i*nrhs+j]=sin(M_PI*(i+1)/(n+1.));
}
}
printf("end ini \n");
int ierr;
// ROW_MAJOR is C order, Lapacke will compute ldab by himself.
ierr=LAPACKE_dgbsv(LAPACK_ROW_MAJOR, n, kl,ku,nrhs, a,n, ipiv, b,nrhs );
if(ierr<0){LAPACKE_xerbla( "LAPACKE_dgbsv", ierr );}
printf("output of LAPACKE_dgbsv\n");
for(i=0;i<n;i++){
for(j=0;j<nrhs;j++){
printf("%g ",b[i*nrhs+j]);
}
printf("\n");
}
//checking correctness
double norm=0;
double diffnorm=0;
for(i=0;i<n;i++){
for(j=0;j<nrhs;j++){
norm+=b[i*nrhs+j]*b[i*nrhs+j];
diffnorm+=(b[i*nrhs+j]-1./(M_PI*M_PI)*sin(M_PI*(i+1)/(n+1.)))*(b[i*nrhs+j]-1./(M_PI*M_PI)*sin(M_PI*(i+1)/(n+1.)));
}
}
printf("analical solution is 1/(PI*PI)*sin(x)\n");
printf("relative difference is %g\n",sqrt(diffnorm/norm));
free(a);
free(b);
free(ipiv);
return 0;
}
我正在尝试使用 LAPACK 的 C 接口求解一般带状矩阵,在 Intel 的 MKL 中称为 LAPACKE。我要调用的函数是 *gbsv
,其中 *
表示格式。不幸的是,我发现 非常 很难找到有关如何使用 C 接口格式化带状矩阵的工作示例。如果有人可以为那里的所有 C 用户提供一个工作示例,我向你保证它会有所帮助。
fortran 布局作为示例给出 here,但我不确定如何将其格式化以输入到 LAPACKE。我还应该注意,在我的问题中,我必须即时构建带状矩阵。所以我有 5 个系数,A、B、C、D、E 用于每个 i 节点,必须将它们放入带状矩阵形式,然后传递给 LAPACKE。
函数LAPACKE_dgbsv()
的原型如下:
lapack_int LAPACKE_dgbsv( int matrix_layout, lapack_int n, lapack_int kl,
lapack_int ku, lapack_int nrhs, double* ab,
lapack_int ldab, lapack_int* ipiv, double* b,
lapack_int ldb )
与Lapack的函数dgbsv()
的主要区别在于参数matrix_layout
,可以是LAPACK_ROW_MAJOR
(C序)或LAPACK_COL_MAJOR
(Fortran序)。如果 LAPACK_ROW_MAJOR
,LAPACKE_dgbsv
将转置矩阵,调用 dgbsv()
然后将矩阵转置回 C 排序。
其他参数的含义与函数dgbsv()
相同。如果使用 LAPACK_ROW_MAJOR
,则 dgbsv()
的正确 ldab
将由 LAPACKE_dgbsv()
计算,参数 ldab
可以设置为 n
。但是,就像 dgbsv()
一样,必须为矩阵 ab
分配额外的 space 来存储因式分解的细节。
下面的例子利用LAPACKE_dgbsv()
通过中心有限差分求解一维稳态扩散。考虑了零温度边界条件,并将正弦波之一用作源项以检查正确性。以下程序由gcc main3.c -o main3 -llapacke -llapack -lblas -Wall
编译:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <lapacke.h>
int main(void){
srand (time(NULL));
//size of the matrix
int n=10;
// number of right-hand size
int nrhs=4;
int ku=2;
int kl=2;
// ldab is larger than the number of bands,
// to store the details of factorization
int ldab = 2*kl+ku+1;
//memory initialization
double *a=malloc(n*ldab*sizeof(double));
if(a==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
double *b=malloc(n*nrhs*sizeof(double));
if(b==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
int *ipiv=malloc(n*sizeof(int));
if(ipiv==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
int i,j;
double fact=1*((n+1.)*(n+1.));
//matrix initialization : the different bands
// are stored in rows kl <= j< 2kl+ku+1
for(i=0;i<n;i++){
a[(0+kl)*n+i]=0;
a[(1+kl)*n+i]=-1*fact;
a[(2+kl)*n+i]=2*fact;
a[(3+kl)*n+i]=-1*fact;
a[(4+kl)*n+i]=0;
//initialize source terms
for(j=0;j<nrhs;j++){
b[i*nrhs+j]=sin(M_PI*(i+1)/(n+1.));
}
}
printf("end ini \n");
int ierr;
// ROW_MAJOR is C order, Lapacke will compute ldab by himself.
ierr=LAPACKE_dgbsv(LAPACK_ROW_MAJOR, n, kl,ku,nrhs, a,n, ipiv, b,nrhs );
if(ierr<0){LAPACKE_xerbla( "LAPACKE_dgbsv", ierr );}
printf("output of LAPACKE_dgbsv\n");
for(i=0;i<n;i++){
for(j=0;j<nrhs;j++){
printf("%g ",b[i*nrhs+j]);
}
printf("\n");
}
//checking correctness
double norm=0;
double diffnorm=0;
for(i=0;i<n;i++){
for(j=0;j<nrhs;j++){
norm+=b[i*nrhs+j]*b[i*nrhs+j];
diffnorm+=(b[i*nrhs+j]-1./(M_PI*M_PI)*sin(M_PI*(i+1)/(n+1.)))*(b[i*nrhs+j]-1./(M_PI*M_PI)*sin(M_PI*(i+1)/(n+1.)));
}
}
printf("analical solution is 1/(PI*PI)*sin(x)\n");
printf("relative difference is %g\n",sqrt(diffnorm/norm));
free(a);
free(b);
free(ipiv);
return 0;
}