Format banded Matrix für LAPACKE

Question

Ich versuche, eine allgemeine Bandmatrix mit der C-Schnittstelle zu LAPACK namens LAPACKE in Intels MKL zu lösen. Die Funktion, die ich aufrufen möchte, ist *gbsv, wobei * das Format bezeichnet. Leider finde ich es SEHR schwierig, Arbeitsbeispiele zu finden, wie man die gebänderte Matrix unter Verwendung der C-Schnittstelle formatiert. Wenn jemand ein funktionierendes Beispiel für alle C-Benutzer da draußen zur Verfügung stellen könnte, versichere ich Ihnen, dass es hilfreich wäre.

Das Fortran-Layout wird als Beispiel here angegeben, aber ich bin mir nicht sicher, wie ich dies für die Eingabe in LAPACKE formatieren würde. Ich sollte auch erwähnen, dass ich bei meinem Problem die Banded Matrix on the fly aufbauen muss. Also habe ich 5 Koeffizienten, A, B, C, D, E für jeden i-Knoten, die in eine gebänderte Matrixform gebracht und dann an LAPACKE übergeben werden müssen.

Answer 1

Der Prototyp der Funktion LAPACKE_dgbsv() ist die folgende:

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) 

Der Hauptunterschied mit der Funktion dgbsv() von Lapack ist das Argument matrix_layout, die LAPACK_ROW_MAJOR (C Ordnung) oder LAPACK_COL_MAJOR (Fortran Ordnung) sein kann. Wenn LAPACK_ROW_MAJOR, LAPACKE_dgbsv die Matrizen transponieren, rufen Sie dgbsv() auf und transponieren Sie die Matrizen dann wieder in C-Reihenfolge. Die Bedeutung der anderen Argumente ist dieselbe wie für die Funktion dgbsv(). Wenn LAPACK_ROW_MAJOR verwendet wird, wird die korrekte ldab für dgbsv() von LAPACKE_dgbsv() berechnet und das Argument ldab kann auf n gesetzt werden. Genau wie dgbsv() muss jedoch zusätzlicher Speicherplatz für die Matrix ab reserviert werden, um die Details der Faktorisierung zu speichern.

Das folgende Beispiel verwendet LAPACKE_dgbsv(), um 1D stationäre Diffusion durch zentrierte endliche Differenz zu lösen. Null-Temperatur-Randbedingungen werden berücksichtigt, und eine Sinuswelle wird als Quellterm verwendet, um die Korrektheit zu überprüfen. Das folgende Programm wird zusammengestellt von 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; 
}