Формат полосовой матрицы для LAPACKE

Question

Я пытаюсь решить общую полосу с матрицей, используя интерфейс C для LAPACK, называемый LAPACKE в Intel MKL. Функция, которую я пытаюсь вызвать, - *gbsv, где * обозначает формат. К сожалению, я нахожу это ОЧЕНЬ трудно найти рабочие примеры о том, как отформатировать матрицу с использованием интерфейса C. Если кто-то может предоставить рабочий пример всем пользователям C, я заверяю вас, что это будет полезно.

Формат fortran приведен в качестве примера here, но я не совсем уверен, как бы отформатировать его для ввода в LAPACKE. Я должен также отметить, что в моей проблеме я должен построить полосатую матрицу на лету. Таким образом, у меня есть 5 коэффициентов: A, B, C, D, E для каждого i-узла, которые должны быть помещены в матричную форму, затем переданы в LAPACKE.

Answer 1

Прототип функции 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 (С упорядоченностью) или LAPACK_COL_MAJOR (Fortran упорядочения). Если LAPACK_ROW_MAJOR, LAPACKE_dgbsv транспонирует матрицы, вызовите dgbsv(), а затем переставьте матрицы обратно на упорядочение C.

Смысл других аргументов тот же, что и для функции dgbsv(). Если используется LAPACK_ROW_MAJOR, то ldab для dgbsv() будет вычислено 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; 
}