2

我正在尝试使用英特尔 MKL 中称为 LAPACKE 的 LAPACK 的 C 接口来解决一般带状矩阵。我试图调用的函数是*gbsv，其中*表示格式。不幸的是，我发现很难找到有关如何使用 C 接口格式化带状矩阵的工作示例。如果有人可以为所有 C 用户提供一个工作示例，我向您保证这会有所帮助。

此处给出了 fortran 布局作为示例，但我不确定如何将其格式化以输入到 LAPACKE。我还应该注意，在我的问题中，我必须即时构建带状矩阵。所以我有 5 个系数，每个 i 节点的 A、B、C、D、E，必须放入带状矩阵形式，然后传递给 LAPACKE。

4

1 回答 1

2

该函数的原型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 )

dgbsv()与Lapack函数的主要区别在于参数matrix_layout，它可以是LAPACK_ROW_MAJOR（C 排序）或LAPACK_COL_MAJOR（Fortran 排序）。如果LAPACK_ROW_MAJOR,LAPACKE_dgbsv将转置矩阵，调用dgbsv()然后将矩阵转置回 C 排序。

其他参数的含义与函数相同dgbsv()。如果LAPACK_ROW_MAJOR使用了，那么正确ldab的 fordgbsv()将被计算LAPACKE_dgbsv()并且参数ldab可以设置为n。然而，就像 一样dgbsv()，必须为矩阵分配额外的空间ab来存储分解的细节。

下面的例子利用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;
}
于 2016-04-02T14:39:46.093 回答