SPSOLVE computes the solution of a sparse system through a sequence of steps: Initialization, ordering to reduce fill-in, symbolic factorization, numeric factorization, and triangular solve. A user code can call individual routines or make use of a one-call interface to perform these steps.
Listed in the table below are user-accessible routines in SPSOLVE and their purposes.
|
Matrices with the same structure but with different numerical values can be solved by calling SPSOLVE routines in the following order shown:
call dgssin() ! initialization, input coefficient matrix structure call dgssor() ! fill-reducing ordering, symbolic factorization ! (or call dgssuo() to specify a user ordering, ! and perform symbolic factorization) do m = 1, number_of_structurally_identical_matrices call dgssfa() ! input coefficient matrix values, numeric ! factorization do r = 1, number_of_right_hand_sides call dgsssl() ! triangular solve enddo enddo
The one-call interface is not as flexible as the regular interface, but it covers the most common case of factoring a single matrix and solving some number of right-hand sides. Additional calls to dgsssl() are used to solve for additional right-hand sides, as shown in the following example.
call dgssfs() ! initialization, input coefficient matrix structure ! fill-reducing ordering, symbolic factorization ! input coefficient matrix values, numeric factorization ! triangular solve do r = 1, number_of_right_hand_sides call dgsssl() ! triangular solve enddo
To use SPSOLVE, you must call its routines in the following order shown:
One-Call Interface: For solving single matrix
DGSSFS() - Initialize, order, factor, solve
DGSSSL() - Additional solves (optional): repeat DGSSSL() as needed
DGSSDA() - Deallocate working storage
Regular Interface: For solving multiple matrices with the same structure
DGSSIN() - Initialize
DGSSOR() or DGSSUO() - Order and symbolically factor
DGSSFA() - Factor
DGSSSL() - Solve: repeat DGSSFA() or DGSSSL() as needed
DGSSDA() - Deallocate working storage
The following examples show solving a symmetric system using the one-call interface, and solving a symmetric system using the regular interface.
In Example 1, the one-call interface is used to solve a symmetric system, and in Example 2, individual routines are called to solve a symmetric system.
Example 5 shows how the Fortran SPSOLVE interface can be called from a C program. For more information on how to call Fortran routines from C programs, see the Oracle Solaris Studio 12.4: Fortran Programming Guide.
Example 1 Solving a Symmetric System-One-Call Interfacemy_system% cat example_1call.f program example_1call c c This program is an example driver that calls the sparse solver. c It factors and solves a symmetric system, by calling the c one-call interface. c implicit none integer neqns, ier, msglvl, outunt, ldrhs, nrhs character mtxtyp*2, pivot*1, ordmthd*3 double precision handle(150) integer colstr(6), rowind(9) double precision values(9), rhs(5), xexpct(5) integer i c c Sparse matrix structure and value arrays. From George and Liu, c page 3. c Ax = b, (solve for x) where: c c 4.0 1.0 2.0 0.5 2.0 2.0 7.0 c 1.0 0.5 0.0 0.0 0.0 2.0 3.0 c A = 2.0 0.0 3.0 0.0 0.0 x = 1.0 b = 7.0 c 0.5 0.0 0.0 0.625 0.0 -8.0 -4.0 c 2.0 0.0 0.0 0.0 16.0 -0.5 -4.0 c data colstr / 1, 6, 7, 8, 9, 10 / data rowind / 1, 2, 3, 4, 5, 2, 3, 4, 5 / data values / 4.0d0, 1.0d0, 2.0d0, 0.5d0, 2.0d0, 0.5d0, 3.0d0, & 0.625d0, 16.0d0 / data rhs / 7.0d0, 3.0d0, 7.0d0, -4.0d0, -4.0d0 / data xexpct / 2.0d0, 2.0d0, 1.0d0, -8.0d0, -0.5d0 / c c set calling parameters c mtxtyp= 'ss' pivot = 'n' neqns = 5 nrhs = 1 ldrhs = 5 outunt = 6 msglvl = 0 ordmthd = 'mmd' c c call single call interface c call dgssfs ( mtxtyp, pivot, neqns , colstr, rowind, & values, nrhs , rhs, ldrhs , ordmthd, & outunt, msglvl, handle, ier ) if ( ier .ne. 0 ) goto 110 c c deallocate sparse solver storage c call dgssda ( handle, ier ) if ( ier .ne. 0 ) goto 110 c c print values of sol c 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 c call to sparse solver returns an error c write ( 6 , 400 ) & ' example: FAILED sparse solver error number = ', ier stop 200 format(a5,3a20) 300 format(i5,3d20.12) ! i/sol/xexpct values 400 format(a60,i20) ! fail message, sparse solver error number end my_system% f95 -dalign example_1call.f -library=sunperf my_sytem% a.out i rhs(i) expected rhs(i) error 1 0.200000000000D+01 0.200000000000D+01 -0.528466159722D-13 2 0.200000000000D+01 0.200000000000D+01 0.105249142734D-12 3 0.100000000000D+01 0.100000000000D+01 0.350830475782D-13 4 -0.800000000000D+01 -0.800000000000D+01 0.426325641456D-13 5 -0.500000000000D+00 -0.500000000000D+00 0.660582699652D-14Example 2 Solving a Symmetric System – Regular Interface
my_system% cat example_ss.f program example_ss c c This program is an example driver that calls the sparse solver. c It factors and solves a symmetric system. implicit none integer neqns, ier, msglvl, outunt, ldrhs, nrhs character mtxtyp*2, pivot*1, ordmthd*3 double precision handle(150) integer colstr(6), rowind(9) double precision values(9), rhs(5), xexpct(5) integer i c c Sparse matrix structure and value arrays. From George and Liu, c page 3. c Ax = b, (solve for x) where: c c 4.0 1.0 2.0 0.5 2.0 2.0 7.0 c 1.0 0.5 0.0 0.0 0.0 2.0 3.0 c A = 2.0 0.0 3.0 0.0 0.0 x = 1.0 b = 7.0 c 0.5 0.0 0.0 0.625 0.0 -8.0 -4.0 c 2.0 0.0 0.0 0.0 16.0 -0.5 -4.0 c data colstr / 1, 6, 7, 8, 9, 10 / data rowind / 1, 2, 3, 4, 5, 2, 3, 4, 5 / data values / 4.0d0, 1.0d0, 2.0d0, 0.5d0, 2.0d0, 0.5d0, & 3.0d0, 0.625d0, 16.0d0 / data rhs / 7.0d0, 3.0d0, 7.0d0, -4.0d0, -4.0d0 / data xexpct / 2.0d0, 2.0d0, 1.0d0, -8.0d0, -0.5d0 / c c initialize solver c mtxtyp= 'ss' pivot = 'n' neqns = 5 outunt = 6 msglvl = 0 c c call regular interface c call dgssin ( mtxtyp, pivot, neqns , colstr, rowind, & outunt, msglvl, handle, ier ) if ( ier .ne. 0 ) goto 110 c c ordering and symbolic factorization c ordmthd = 'mmd' call dgssor ( ordmthd, handle, ier ) if ( ier .ne. 0 ) goto 110 c c numeric factorization c call dgssfa ( neqns, colstr, rowind, values, handle, ier ) if ( ier .ne. 0 ) goto 110 c c solution c nrhs = 1 ldrhs = 5 call dgsssl ( nrhs, rhs, ldrhs, handle, ier ) if ( ier .ne. 0 ) goto 110 c c deallocate sparse solver storage c call dgssda ( handle, ier ) if ( ier .ne. 0 ) goto 110 c c print values of sol c 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 c call to sparse solver returns an error c write ( 6 , 400 ) & ' example: FAILED sparse solver error number = ', ier stop 200 format(a5,3a20) 300 format(i5,3d20.12) ! i/sol/xexpct values 400 format(a60,i20) ! fail message, sparse solver error number end my_system% f95 -dalign example_ss.f -library=sunperf my_sytem% a.out i rhs(i) expected rhs(i) error 1 0.200000000000D+01 0.200000000000D+01 -0.528466159722D-13 2 0.200000000000D+01 0.200000000000D+01 0.105249142734D-12 3 0.100000000000D+01 0.100000000000D+01 0.350830475782D-13 4 -0.800000000000D+01 -0.800000000000D+01 0.426325641456D-13 5 -0.500000000000D+00 -0.500000000000D+00 0.660582699652D-14Example 3 Solving a Structurally Symmetric System With Unsymmetric Values – Regular Interface
my_system% cat example_su.f program example_su c c This program is an example driver that calls the sparse solver. c It factors and solves a structurally symmetric system c (w/unsymmetric values). c implicit none integer neqns, ier, msglvl, outunt, ldrhs, nrhs character mtxtyp*2, pivot*1, ordmthd*3 double precision handle(150) integer colstr(5), rowind(8) double precision values(8), rhs(4), xexpct(4) integer i c c Sparse matrix structure and value arrays. Coefficient matrix c has a symmetric structure and unsymmetric values. c Ax = b, (solve for x) where: c c 1.0 3.0 0.0 0.0 1.0 7.0 c 2.0 4.0 0.0 7.0 2.0 38.0 c A = 0.0 0.0 6.0 0.0 x = 3.0 b = 18.0 c 0.0 5.0 0.0 8.0 4.0 42.0 c data colstr / 1, 3, 6, 7, 9 / data rowind / 1, 2, 1, 2, 4, 3, 2, 4 / data values / 1.0d0, 2.0d0, 3.0d0, 4.0d0, 5.0d0, 6.0d0, 7.0d0, & 8.0d0 / data rhs / 7.0d0, 38.0d0, 18.0d0, 42.0d0 / data xexpct / 1.0d0, 2.0d0, 3.0d0, 4.0d0 / c c initialize solver c mtxtyp= 'su' pivot = 'n' neqns = 4 outunt = 6 msglvl = 0 c c call regular interface c call dgssin ( mtxtyp, pivot, neqns , colstr, rowind, & outunt, msglvl, handle, ier ) if ( ier .ne. 0 ) goto 110 c c ordering and symbolic factorization c ordmthd = 'mmd' call dgssor ( ordmthd, handle, ier ) if ( ier .ne. 0 ) goto 110 c c numeric factorization c call dgssfa ( neqns, colstr, rowind, values, handle, ier ) if ( ier .ne. 0 ) goto 110 c c solution c nrhs = 1 ldrhs = 4 call dgsssl ( nrhs, rhs, ldrhs, handle, ier ) if ( ier .ne. 0 ) goto 110 c c deallocate sparse solver storage c call dgssda ( handle, ier ) if ( ier .ne. 0 ) goto 110 c c print values of sol c 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 c call to sparse solver returns an error c write ( 6 , 400 ) & ' example: FAILED sparse solver error number = ', ier stop 200 format(a5,3a20) 300 format(i5,3d20.12) ! i/sol/xexpct values 400 format(a60,i20) ! fail message, sparse solver error number end my_system% f95 -dalign example_su.f -library=sunperf my_system% a.out i rhs(i) expected rhs(i) error 1 0.100000000000D+01 0.100000000000D+01 0.000000000000D+00 2 0.200000000000D+01 0.200000000000D+01 0.000000000000D+00 3 0.300000000000D+01 0.300000000000D+01 0.000000000000D+00 4 0.400000000000D+01 0.400000000000D+01 0.000000000000D+00Example 4 Solving an Unsymmetric System – Regular Interface
my_system% cat example_uu.f program example_uu c c This program is an example driver that calls the sparse solver. c It factors and solves an unsymmetric system. c implicit none integer neqns, ier, msglvl, outunt, ldrhs, nrhs character mtxtyp*2, pivot*1, ordmthd*3 double precision handle(150) integer colstr(6), rowind(10) double precision values(10), rhs(5), xexpct(5) integer i c c Sparse matrix structure and value arrays. Unsummetric matrix A. c Ax = b, (solve for x) where: c c 1.0 0.0 0.0 0.0 0.0 1.0 1.0 c 2.0 6.0 0.0 0.0 9.0 2.0 59.0 c A = 3.0 0.0 7.0 0.0 0.0 x = 3.0 b = 24.0 c 4.0 0.0 0.0 8.0 0.0 4.0 36.0 c 5.0 0.0 0.0 0.0 10.0 5.0 55.0 c data colstr / 1, 6, 7, 8, 9, 11 / data rowind / 1, 2, 3, 4, 5, 2, 3, 4, 2, 5 / data values / 1.0d0, 2.0d0, 3.0d0, 4.0d0, 5.0d0, 6.0d0, 7.0d0, & 8.0d0, 9.0d0, 10.0d0 / data rhs / 1.0d0, 59.0d0, 24.0d0, 36.0d0, 55.0d0 / data xexpct / 1.0d0, 2.0d0, 3.0d0, 4.0d0, 5.0d0 / c c initialize solver c mtxtyp= 'uu' pivot = 'n' neqns = 5 outunt = 6 msglvl = 3 call dgssin ( mtxtyp, pivot, neqns , colstr, rowind, & outunt, msglvl, handle, ier ) if ( ier .ne. 0 ) goto 110 c c ordering and symbolic factorization c ordmthd = 'mmd' call dgssor ( ordmthd, handle, ier ) if ( ier .ne. 0 ) goto 110 c c numeric factorization c call dgssfa ( neqns, colstr, rowind, values, handle, ier ) if ( ier .ne. 0 ) goto 110 c c solution c nrhs = 1 ldrhs = 5 call dgsssl ( nrhs, rhs, ldrhs, handle, ier ) if ( ier .ne. 0 ) goto 110 c c deallocate sparse solver storage c call dgssda ( handle, ier ) if ( ier .ne. 0 ) goto 110 c c print values of sol c 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 c call to sparse solver returns an error c write ( 6 , 400 ) & ' example: FAILED sparse solver error number = ', ier stop 200 format(a5,3a20) 300 format(i5,3d20.12) ! i/sol/xexpct values 400 format(a60,i20) ! fail message, sparse solver error number end my_system% f95 -dalign example_uu.f -library=sunperf my_system% a.out i rhs(i) expected rhs(i) error 1 0.100000000000D+01 0.100000000000D+01 0.000000000000D+00 2 0.200000000000D+01 0.200000000000D+01 0.000000000000D+00 3 0.300000000000D+01 0.300000000000D+01 0.000000000000D+00 4 0.400000000000D+01 0.400000000000D+01 0.000000000000D+00 5 0.500000000000D+01 0.500000000000D+01 0.000000000000D+00Example 5 Calling SPSOLVE Routines from C
#include <stdio.h> #include <stdlib.h> #include <math.h> #include <sys/time.h> #include <sunperf.h> int main() { /* Sparse matrix structure and value arrays. Coefficient matrix is a general unsymmetric sparse matrix. Ax = b, (solve for x) where: 1.0 0.0 7.0 9.0 0.0 1.0 17.0 2.0 4.0 0.0 0.0 0.0 1.0 6.0 A = 0.0 5.0 8.0 0.0 0.0 x = 1.0 b = 13.0 0.0 0.0 0.0 10.0 11.0 1.0 21.0 3.0 6.0 0.0 0.0 12.0 1.0 21.0 */ /* Array indices must be one-based for calling SPSOLVE routines */ int colstr[] = {1, 4, 7, 9, 11, 13}; int rowind[] = {1, 2, 5, 2, 3, 5, 1, 3, 1, 4, 4, 5}; double 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}; double rhs[] = {17.0, 6.0, 13.0, 21.0, 21.0}; double xexpct[] = {1.0, 1.0, 1.0, 1.0, 1.0}; int n = 5, nnz = 12, nrhs = 1, msglvl = 0, outunt = 6, ierr, i,j,k, int_ierr; double t[4], handle[150]; char type[] = "uu", piv = 'n'; /* Last two parameters in argument list indicate lengths of * character arguments type and piv */ dgssin_(type, &piv, &n, colstr, rowind, &outunt, &msglvl, handle, &ierr,2,1); if (ierr != 0) { int_ierr = ierr; printf("dgssin err = %d\n", int_ierr); return -1; } char ordmth[] = "mmd"; dgssor_(ordmth, handle, &ierr, 3); if (ierr != 0) { int_ierr = ierr; printf("dgssor err = %d\n", int_ierr); return -1; } dgssfa_(&n, colstr, rowind, values, handle, &ierr); if (ierr != 0) { int_ierr = ierr; printf("dgssfa err = %d\n", int_ierr); return -1; } dgsssl_(&nrhs, rhs, &n, handle, &ierr); if (ierr != 0) { int_ierr = ierr; printf("dgsssl err = %d\n", int_ierr); return -1; } printf("i computed solution expected solution\n"); for (i=0; i<n; i++) printf("%d %lf %lf\n", i,rhs[i], 1.0); } my_system% cc -m32 -xmemalign=8s dr.c -library=sunperf my_system% ./a.out i computed solution expected solution 0 1.000000 1.000000 1 1.000000 1.000000 2 1.000000 1.000000 3 1.000000 1.000000 4 1.000000 1.000000