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-14
Example 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-14
Example 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+00
Example 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+00
Example 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