LAPACKE

Question

için bantlanmış matris C'nin Intel'in MKL'sindeki LAPACKE adlı LAPACK için C arabirimini kullanarak genel bantlı matrisi çözmeye çalışıyorum. Aramaya çalıştığım işlev, * biçimini belirtildiği *gbsv'dur. Ne yazık ki, C arayüzünü kullanarak bantlı matrisin nasıl biçimlendirileceğine dair çalışma örnekleri bulmak VERY buluyorum. Birisi orada tüm C kullanıcıları için bir çalışma örneği sağlayabilirse, size yardımcı olacağını garanti ederim.

Fortran düzeni, örnek olarak here olarak verilir, ancak LAPACKE giriş için bunu nasıl biçimlendireceğimi tam olarak bilmiyorum. Dikkat etmeliyim ki, benim sorunumda, bantlanmış matrisi anında oluşturmak zorundayım. Bu yüzden, her bir i-düğümü için 5, A, B, C, D, E, yani bir bantlanmış matris formuna konması gereken ve sonra LAPACKE'e geçirilen 5 katsayıya sahibim.

Answer 1

fonksiyonu LAPACKE_dgbsv() prototipidir aşağıdaki gibidir:

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) 

fonksiyonu LAPACK arasında dgbsv() ile temel fark LAPACK_ROW_MAJOR (Cı sipariş) ya da LAPACK_COL_MAJOR (Fortran sıralama) olabilir argüman matrix_layout vardır. LAPACK_ROW_MAJOR, LAPACKE_dgbsv, matrisleri aktarırsa, dgbsv() numaralı telefonu arayın ve ardından matrisleri C sıralamasına geri aktarın.

Diğer bağımsız değişkenlerin anlamı, dgbsv() işleviyle aynıdır. LAPACK_ROW_MAJOR kullanılıyorsa, dgbsv() için doğru ldabLAPACKE_dgbsv() hesaplanır ve ldab argümanı n olarak ayarlanabilir. Ancak, dgbsv() gibi, faktörizasyonun ayrıntılarını saklamak için ab matrisi için ek alan ayrılmalıdır.

Aşağıdaki örnek, merkezileştirilmiş sonlu farkla 1D istasyon difüzyonunu çözmek için LAPACKE_dgbsv()'u kullanır. Sıfır sıcaklık sınır koşulu dikkate alınır ve doğruluk kontrol etmek için bir sinüs dalgasından biri kaynak olarak kullanılır. Aşağıdaki programı gcc main3.c -o main3 -llapacke -llapack -lblas -Wall tarafından derlenen:

benim tesisatın mkl kök klasöründe

#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; 
}