SuperLU インタフェース

言語:

SuperLU には 2 つのドライバルーチン (単純およびエキスパート) があります。これらを呼び出して、SPSOLVE のワンコールインタフェースと同様の方法で、一般非対称スパース連立方程式を完全に解くことができます。これらおよびほかの SuperLU のユーザー呼び出し可能ルーチンは、単精度、倍精度、複素数、および倍精度複素数のデータ型で提供されています。すべての外部ルーチンの単精度名を次の表に示します。これらのルーチンにはマニュアルページ (セクション 3P) が用意されています。また、アプリケーションで使用されているスパース行列データ構造については、SuperMatrix(3P) のマニュアルページを参照してください。

表 5 SuperLU 計算ルーチン

ルーチン	説明
`sgstrf`	因子分解を計算します
`sgssvx`	因子分解と求解 (エキスパートドライバ)
`sgssv`	因子分解と求解 (単純ドライバ)
`sgstrs`	三角求解を計算します
`sgsrfs`	計算解を改良し、誤差限界を提供します
`slangs`	1 ノルム、フロベニウスノルム、または無限大ノルムを計算します
`sgsequ`	行および列のスケーリングを計算します
`sgscon`	条件数の逆数を推定します
`slaqgs`	一般スパース行列を均衡化します

表 6 SuperLU ユーティリティールーチン

ルーチン	説明
`LUSolveTime`	求解段階で費やされた時間を返します
`LUFactTime`	因子分解段階で費やされた時間を返します
`LUFactFlops`	因子分解段階での浮動小数点演算の回数を返します
`LUSolveFlops`	求解段階での浮動小数点演算の回数を返します
`sQuerySpace`	メモリー統計の情報を返します
`sp_ienv`	指定されたマシンに依存するパラメータを返します
`sPrintPerf`	計算ルーチンによって収集された統計を出力します
`set_default_options`	ソルバーの動作を制御するパラメータをデフォルトオプションに設定します
`StatInit`	パフォーマンス統計を格納する構造体を割り当て、初期化します
`StatFree`	パフォーマンス統計を格納する構造体を解放します
`Destroy_Dense_Matrix`	密形式の SuperMatrix を割り当て解除します
`Destroy_SuperNode_Matrix`	スーパーノード形式の SuperMatrix を割り当て解除します
`Destroy_CompCol_Matrix`	圧縮スパース列形式の SuperMatrix を割り当て解除します
`Destroy_CompCol_Permuted`	置換された圧縮スパース列形式の SuperMatrix を割り当て解除します
`Destroy_SuperMatrix_Store`	SuperMatrix の行列を格納する実際の格納領域を割り当て解除します
`sCopy_CompCol_Matrix`	圧縮スパース列形式の SuperMatrix をコピーします
`sCreate_CompCol_Matrix`	圧縮スパース列形式の SuperMatrix を割り当てます
`sCreate_Dense_Matrix`	密形式の SuperMatrix を割り当てます
`sCreate_CompRow_Matrix`	圧縮スパース行形式の SuperMatrix を割り当てます
`sCreate_SuperNode_Matrix`	スーパーノード形式の SuperMatrix を割り当てます
`sp_preorder`	元のスパース行列の列を置換します
`sp_sgemm`	SuperMatrix に密行列を乗算します

C からの SuperLU の呼び出し

SuperLU のルーチンは C で記述されています。したがって、列および行に関連するインデックスは 0 オリジンでなければなりません。次の例では、倍精度の単純ドライバ dgssv() を呼び出して、因子 L および U を計算し、解行列を求めます。

使用例 6 SuperLU 単純ドライバ

#include <stdio.h>
#include <sunperf.h>
 
#define M 5
#define N 5
 
int main(int argc, char *argv[])
{
    SuperMatrix A, L, U, B1, B2;
    int      perm_r[M];  /* row permutations from partial pivoting */
    int      perm_c[N];  /* column permutation vector */
    int      info, i;
    superlu_options_t options;
    SuperLUStat_t stat;
    trans_t  trans = NOTRANS;
 
    printf("Example code calling SuperLU simple driver to factor a \n");
    printf("general unsymmetric matrix and solve two right-hand-side matrices\n");
 
   /* the matrix in Harwell-Boeing format. */
    int m = M;
    int n = M;
    int nnz = 12;
    double *dp;
   /* nonzeros of A, column-wise */
    double a[] = {1.0, 2.0, 3.0,  4.0, 5.0,  6.0,
                  7.0, 8.0, 9.0, 10.0, 11.0, 12.0};
   /* row index of nonzeros */
    int asub[] = {0, 1, 4, 1, 2, 4, 0, 2, 0, 3, 3, 4};
   /* column pointers */
    int xa[]   = {0, 3, 6, 8, 10, 12};
 
    /* Create Matrix A in the format expected by SuperLU */
    dCreate_CompCol_Matrix(&A, m, n, nnz, a, asub, xa, SLU_NC, SLU_D, SLU_GE);
 
    int nrhs = 1;
    double rhs1[] = {17.0, 6.0, 13.0, 21.0, 21.0};
    double rhs2[] = {17*.3, 6*.3, 13*.3, 21*.3, 21*.3};
 
    /* right-hand side matrix B1, B2 */
    dCreate_Dense_Matrix(&B1, m, nrhs, rhs1, m, SLU_DN, SLU_D, SLU_GE);
    dCreate_Dense_Matrix(&B2, m, nrhs, rhs2, m, SLU_DN, SLU_D, SLU_GE);
 
    /* set options that control behavior of solver to default parameters */
    set_default_options(&options);
    options.ColPerm = NATURAL;
 
    /* Initialize the statistics variables. */
    StatInit(&stat);
    /* factor input matrix and solve the first right-hand-side matrix */
    dgssv(&options, &A, perm_c, perm_r, &L, &U, &B1, &stat, &info);
 
    printf("\nsolution matrix B1:\n");
    dp = (double *) (((NCformat *)B1.Store)->nzval);
    printf("    i    rhs[i]     expected\n");
    for (i=0; i<M; i++)
      printf("%5d   %7.4lf     %7.4lf\n", i, dp[i], 1.0);
    printf("Factor time  = %8.2e sec\n", stat.utime[FACT]);
    printf("Solve time   = %8.2e sec\n\n\n", stat.utime[SOLVE]);
 
    /* solve the second right-hand-side matrix */
    dgstrs(trans, &L, &U, perm_c, perm_r, &B2, &stat, &info);
 
    printf("solution matrix B2:\n");
    dp = (double *) (((NCformat *)B2.Store)->nzval);
    printf("    i    rhs[i]     expected\n");
    for (i=0; i<M; i++)
      printf("%5d   %7.4lf     %7.4lf\n", i, dp[i], 0.3);
    printf("Solve time   = %8.2e sec\n", stat.utime[SOLVE]);
 
    StatFree(&stat);
    Destroy_CompCol_Matrix(&A);
    Destroy_SuperMatrix_Store(&B1);
    Destroy_SuperMatrix_Store(&B2);
    Destroy_SuperNode_Matrix(&L);
    Destroy_CompCol_Matrix(&U);
}

上記の例を実行すると、次のようになります。

my_system% cc -xmemalign=8s simple.c -library=sunperf
my_system% a.out
 
Example code calling SuperLU simple driver to factor a
general unsymmetric matrix and solve two right-hand-side matrices
 
solution matrix B1:
    i    rhs[i]     expected
    0    1.0000      1.0000
    1    1.0000      1.0000
    2    1.0000      1.0000
    3    1.0000      1.0000
    4    1.0000      1.0000
Factor time  = 5.43e-02 sec
Solve time   = 6.76e-03 sec
 
 
solution matrix B2:
    i    rhs[i]     expected
    0    0.3000      0.3000
    1    0.3000      0.3000
    2    0.3000      0.3000
    3    0.3000      0.3000
    4    0.3000      0.3000
Solve time   = 6.76e-03 sec

使用例 7 SuperLU エキスパートドライバ

#include <stdio.h>
#include <sunperf.h>
 
#define M    5
#define N    5
#define NRHS 1
 
int main(int argc, char *argv[])
{
    SuperMatrix A, L, U, B, X;
    int      perm_r[M];  /* row permutations from partial pivoting */
    int      perm_c[N];  /* column permutation vector */
    int      etree[N];   /* elimination tree */
    double   ferr[NRHS]; /* estimated forward error bound */
    double   berr[NRHS]; /* component-wise relative backward error */
    double   C[N], R[M]; /* column and row scale factors */
    double   rpg, rcond;
    char     equed[1];   /* Specifies the form of equilibration that was done */          
    double   *work, *dp; /* user-supplied workspace */
    int      lwork = 0;  /* 0 for workspace to be allocated by system malloc */
    int      info, i;
    superlu_options_t options;
    SuperLUStat_t stat;
    mem_usage_t    mem_usage;
 
    printf("Example code calling SuperLU expert driver\n\n");
 
    /* the matrix in Harwell-Boeing format. */
    int m = M;
    int n = M;
    int nnz = 12;
   /* nonzeros of A, column-wise */
    double a[] = {1.0, 2.0, 3.0,  4.0, 5.0,  6.0,
                  7.0, 8.0, 9.0, 10.0, 11.0, 12.0};
    /* row index of nonzeros */
    int asub[] = {0, 1, 4, 1, 2, 4, 0, 2, 0, 3, 3, 4};
    /* column pointers */
    int xa[]   = {0, 3, 6, 8, 10, 12};
    int nrhs = NRHS;
    double rhs[] = {17.0, 6.0, 13.0, 21.0, 21.0};
 
    /* Create Matrix A in the format expected by SuperLU */
    dCreate_CompCol_Matrix(&A, m, n, nnz, a, asub, xa, SLU_NC, SLU_D, SLU_GE);
 
    /* right-hand-side matrix B */
    dCreate_Dense_Matrix(&B, m, nrhs, rhs, m, SLU_DN, SLU_D, SLU_GE);
 
    /*  solution matrix X */
    dCreate_Dense_Matrix(&X, m, nrhs, rhs, m, SLU_DN, SLU_D, SLU_GE);
    set_default_options(&options);
    options.ColPerm = NATURAL;
 
    /* Initialize the statistics variables. */
    StatInit(&stat);
 
    dgssvx(&options, &A, perm_c, perm_r, etree, equed, R, C, &L, &U, work, lwork,
           &B, &X, &rpg, &rcond, ferr, berr, &mem_usage, &stat, &info);
    dp = (double *) (((NCformat *)X.Store)->nzval);
    printf("    i    rhs[i]     expected\n");
    for (i=0; i<M; i++)
      printf("%5d   %7.4lf     %7.4lf\n",
             i, dp[i], 1.0);
    printf("Factor time  = %8.2e sec\n", stat.utime[FACT]);
    printf("Solve time   = %8.2e sec\n", stat.utime[SOLVE]);
 
    StatFree(&stat);
    Destroy_CompCol_Matrix(&A);
    Destroy_SuperMatrix_Store(&B);
    Destroy_SuperNode_Matrix(&L);
    Destroy_CompCol_Matrix(&U);
}

上記の例を実行すると、次のようになります。

my_system% cc -xmemalign=8s expert.c -library=sunperf
my_system% a.out
Example code calling SuperLU expert driver
 
    i    rhs[i]     expected
    0    1.0000      1.0000
    1    1.0000      1.0000
    2    1.0000      1.0000
    3    1.0000      1.0000
    4    1.0000      1.0000
Factor time  = 1.25e-03 sec
Solve time   = 1.70e-04 sec

Fortran からの SuperLU の呼び出し

Fortran から SuperLU を呼び出すもっとも簡単な方法は、SPSOLVE インタフェースを使用することです。SuperLU は、非対称係数行列を解くためにルーチン DGSSIN() の入力引数 MTXTYP を介して選択できます。これは SPSOLVE の初期化ルーチンです。同じ引数がワンコールインタフェースルーチン DGSSFS() にも存在します。

MTXTYP の有効なオプションを次の表に示します。SuperLU を呼び出すには、行列タイプとして 's0' または 'S0' を選択します。SPSOLVE は Fortran ベースなので、入力行列に関連する列および行のインデックスはすべて 1 オリジンであるべきです。ただし、SuperLU を DGSSIN() または DGSSFS() から呼び出す (MTXTYP = 's0' または 'S0' を設定) 場合、これらのインデックスは 0 オリジンでなければなりません。

表 7 DGSSIN() および DGSSFS() の行列タイプオプション

オプション	行列のタイプ	ソルバー
'sp' または 'SP'	対称構造、正定値	SPSOLVE
'ss' または 'SS'	対称構造、対称値	SPSOLVE
'su' または 'SU'	対称構造、非対称値	SPSOLVE
'uu' または 'UU'	非対称構造、非対称値	SPSOLVE
's0' または 'S0'	非対称構造、非対称値	SuperLU

DGSSIN() のあとにルーチン DGSSOR() を呼び出して、フィルインを低減する順序付けおよび記号分解を実行する必要があります。文字引数 (ORDMTHD) は、順序付け方法を選択するために使用されます。この引数はワンコールインタフェースルーチン DGSSFS() にも存在します。SPSOLVE および SuperLU の有効な順序付け方法を次の表に示します。DGSSOR() の代わりに DGSSUO() を呼び出すことにより、特定の順序付けをソルバーに提供することもできます。入力の置換配列は 0 オリジンでなければなりません。

表 8 DGSSOR() および DGSSFS() の行列順序付けオプション

オプション	順序付け方法	ソルバー
'nat' または 'NAT'	自然順序 (順序付けなし)	SPSOLVE、SuperLU
'mmd' または 'MMD'	A'*A の最小次数 (デフォルト)	SPSOLVE、SuperLU
'gnd' または 'GND'	一般入れ子分割	SPSOLVE
'spm' または 'SPM'	A'+A の最小次数順序付け	SuperLU
'sam' または 'SAM'	近似最小次数列	SuperLU

上記のとおり、一般入れ子分割法は SuperLU では使用できません。その一方で、A'+A の最小次数順序付けおよび近似最小次数列順序付けは SPSOLVE では使用できません。

SuperLU の例

次のコード例では、SPSOLVE の通常インタフェースおよびワンコールインタフェースを介して SuperLU を選択し、一般非対称行列の連立方程式の因子分解および求解を行う方法を示します。

使用例 8 SPSOLVE の通常インタフェースを介して SuperLU を呼び出す

      program SLU
 
c  This program is an example driver that calls the regular interface of SPSOLVE
c  to invoke SuperLU to factor and solve a general unsymmetric system.
 
      implicit none
      integer           neqns, ier, msglvl, outunt, ldrhs, nrhs, i
      character         mtxtyp*2, pivot*1, ordmthd*3
      double precision  handle(150)
      integer           colstr(6), rowind(12)
      double precision  values(12), rhs(5), xexpct(5)
 
c  Sparse matrix structure and value arrays.  Coefficient matrix
c    is a general unsymmetric sparse matrix.
c    Ax = b, (solve for x) where:
 
c       1.0   0.0   7.0   9.0   0.0      1.0       17.0
c       2.0   4.0   0.0   0.0   0.0      1.0        6.0
c  A =  0.0   5.0   8.0   0.0   0.0  x = 1.0   b = 13.0
c       0.0   0.0   0.0  10.0  11.0      1.0       21.0
c       3.0   6.0   0.0   0.0  12.0      1.0       21.0
 
c  Array indices must be zero-based for calling SuperLU
      data colstr / 0, 3, 6, 8, 10, 12 /
      data rowind / 0, 1, 4, 1, 2, 4, 0, 2, 0, 3, 3, 4 /
      data values / 1.0, 2.0, 3.0,  4.0, 5.0,  6.0,
     $              7.0, 8.0, 9.0, 10.0, 11.0, 12.0 /
      data rhs    / 17.0, 6.0, 13.0, 21.0, 21.0 /
      data xexpct / 1.0d0, 1.0d0, 1.0d0, 1.0d0, 1.0d0 /
 
c  initialize solver
      mtxtyp= 's0'
      pivot = 'n'
      neqns  = 5
      outunt = 6
      msglvl = 0
 
c  call regular interface
      call dgssin(mtxtyp, pivot,  neqns , colstr, rowind,outunt, msglvl,
     &            handle, ier)
      if ( ier .ne. 0 ) goto 110
 
c  ordering and symbolic factorization
      ordmthd = 'mmd'
      call dgssor(ordmthd, handle, ier)
      if ( ier .ne. 0 ) goto 110
 
c  numeric factorization
      call dgssfa ( neqns, colstr, rowind, values, handle, ier )
      if ( ier .ne. 0 ) goto 110
 
c  solution
      nrhs   = 1
      ldrhs  = 5
      call dgsssl ( nrhs, rhs, ldrhs, handle, ier )
      if ( ier .ne. 0 ) goto 110
 
c  deallocate sparse solver storage
      call dgssda ( handle, ier )
      if ( ier .ne. 0 ) goto 110
 
c  print values of sol
      write(6,200) 'i', 'rhs(i)', 'expected rhs(i)', 'error'
      do i = 1, neqns
        write(6,300) i, rhs(i), xexpct(i), (rhs(i)-xexpct(i))
      enddo
      stop
 
  110 continue
c call to sparse solver returns an error
      write ( 6 , 400 )
     &      ' example: FAILED sparse solver error number = ', ier
      stop
 
  200 format(4x,a1,3x,a6,3x,a15,4x,a6)
  300 format(i5,3x,f5.2,7x,f5.2,8x,e10.2)     ! i/sol/xexpct values
  400 format(a60,i20)   ! fail message, sparse solver error number
      end

上記の例を実行すると、次のようになります。

my_system% f95 -dalign slu.f -library=sunperf
my_system% a.out
 
    i   rhs(i)   expected rhs(i)     error
    1    1.00        1.00          0.00E+00
    2    1.00        1.00         -0.33E-15
    3    1.00        1.00          0.22E-15
    4    1.00        1.00         -0.11E-15
    5    1.00        1.00          0.22E-15

使用例 9 SPSOLVE のワンコールインタフェースを介して SuperLU を呼び出す

program SLU_SINGLE
c  This program is an example driver that calls the regular interface of SPSOLVE
c  to invoke SuperLU to factor and solve a general unsymmetric system.
 
      implicit none
      integer           neqns, ier, msglvl, outunt, ldrhs, nrhs, i
      character         mtxtyp*2, pivot*1, ordmthd*3
      double precision  handle(150)
      integer           colstr(6), rowind(12)
      double precision  values(12), rhs(5), xexpct(5)
 
c  Sparse matrix structure and value arrays.  Coefficient matrix
c    is a general unsymmetric sparse matrix.
c    Ax = b, (solve for x) where:
 
c       1.0   0.0   7.0   9.0   0.0      1.0       17.0
c       2.0   4.0   0.0   0.0   0.0      1.0        6.0
c  A =  0.0   5.0   8.0   0.0   0.0  x = 1.0   b = 13.0
c       0.0   0.0   0.0  10.0  11.0      1.0       21.0
c       3.0   6.0   0.0   0.0  12.0      1.0       21.0
 
c  Array indices must be zero-based for calling SuperLU
      data colstr / 0, 3, 6, 8, 10, 12 /
      data rowind / 0, 1, 4, 1, 2, 4, 0, 2, 0, 3, 3, 4 /
      data values / 1.0, 2.0, 3.0,  4.0, 5.0,  6.0,
     $              7.0, 8.0, 9.0, 10.0, 11.0, 12.0 /
      data rhs    / 17.0, 6.0, 13.0, 21.0, 21.0 /
      data xexpct / 1.0d0, 1.0d0, 1.0d0, 1.0d0, 1.0d0 /
 
c  initialize solver
      mtxtyp= 's0'
      pivot = 'n'
      neqns  = 5
      outunt = 6
      msglvl = 0
      ordmthd = 'mmd'
      nrhs = 1
      ldrhs = 5
 
c  One-call routine of SPSOLVE
      call dgssfs (mtxtyp, pivot, neqns , colstr, rowind,
     &             values, nrhs , rhs, ldrhs , ordmthd,
     &             outunt, msglvl, handle, ier)
      if ( ier .ne. 0 ) goto 110
 
c  deallocate sparse solver storage
      call dgssda ( handle, ier )
      if ( ier .ne. 0 ) goto 110
 
c  print values of sol
      write(6,200) 'i', 'rhs(i)', 'expected rhs(i)', 'error'
      do i = 1, neqns
        write(6,300) i, rhs(i), xexpct(i), (rhs(i)-xexpct(i))
      enddo
      stop
 
  110 continue
c call to sparse solver returns an error
      write ( 6 , 400 )
     &      ' example: FAILED sparse solver error number = ', ier
      stop

  200 format(4x,a1,3x,a6,3x,a15,4x,a6)
  300 format(i5,3x,f5.2,7x,f5.2,8x,e10.2)     ! i/sol/xexpct values
  400 format(a60,i20)   ! fail message, sparse solver error number
      end

上記の例を実行すると、次のようになります。

my_system% f95 -dalign slu_single.f -library=sunperf
my_system% a.out

    i   rhs(i)   expected rhs(i)     error
    1    1.00        1.00          0.00E+00
    2    1.00        1.00         -0.33E-15
    3    1.00        1.00          0.22E-15
    4    1.00        1.00         -0.11E-15
    5    1.00        1.00          0.22E-15