SuperLU has two driver routines, simple and expert, that can be called to completely solve a general unsymmetric sparse system in a similar manner to the one-call interface in SPSOLVE. These and other SuperLU user-callable routines are available in single precision, double precision, complex and double complex data types. Single precision names of all external routines are listed in the following tables. Man pages (section 3P) are available for these routines. Also see the man page of SuperMatrix(3P) for a description of the sparse matrix data structure that is used in the application.
|
|
SuperLU routines are written in C. Therefore, column- and row-related indices must be zero-based. In the following example, double precision simple driver dgssv() is called to compute factors L and U and to solve for the solution matrix.
Example 6 SuperLU Simple Driver#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);
}
Running the above example:
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
Example 7 SuperLU Expert Driver
#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);
}
Running the above example:
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
The simplest way to call SuperLU from Fortran is through the SPSOLVE interface. SuperLU can be selected to solve an unsymmetric coefficient matrix through input argument MTXTYP of routine DGSSIN(), which is the initialization routine in SPSOLVE. The same argument also exists in the one-call interface routine DGSSFS().
Valid options for MTXTYP are listed in the following table. To invoke SuperLU, select 's0' or 'S0' as matrix type. Since SPSOLVE is Fortran-based, all column and row indices associated with the input matrix should be one-based. However, if SuperLU is invoked through DGSSIN() or DGSSFS() (by setting MTXTYP = 's0' or 'S0'), these indices must be zero-based.
|
A call to routine DGSSOR() must follow DGSSIN() to perform fill-reduce ordering and symbolic factorization. A character argument (ORDMTHD) is used to select the desired ordering method. This argument also exists in the one-call interface routine DGSSFS(). Valid ordering methods for SPSOLVE and SuperLU are listed in the following table. You can also provide a particular ordering to the solver by calling DGSSUO() in place of DGSSOR(). The input permutation array must be zero-based.
|
As shown above, the general nested dissection method is not available in SuperLU. On the other hand, the minimum degree ordering on A'+A and approximate minimum degree column ordering are not available in SPSOLVE.
The following code examples show how SuperLU can be selected through the regular interface and the one-call interface of SPSOLVE to factorize and solve a general unsymmetric system of equations.
Example 8 Invoking SuperLU Through SPSOLVE Regular Interface 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
Running the above example:
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
Example 9 Invoking SuperLU through One-Call SPSOLVE Interface
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
Running the above example:
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