00001 /* dspsv.f -- translated by f2c (version 20061008). 00002 You must link the resulting object file with libf2c: 00003 on Microsoft Windows system, link with libf2c.lib; 00004 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 00005 or, if you install libf2c.a in a standard place, with -lf2c -lm 00006 -- in that order, at the end of the command line, as in 00007 cc *.o -lf2c -lm 00008 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 00009 00010 http://www.netlib.org/f2c/libf2c.zip 00011 */ 00012 00013 #include "f2c.h" 00014 #include "blaswrap.h" 00015 00016 /* Subroutine */ int dspsv_(char *uplo, integer *n, integer *nrhs, doublereal 00017 *ap, integer *ipiv, doublereal *b, integer *ldb, integer *info) 00018 { 00019 /* System generated locals */ 00020 integer b_dim1, b_offset, i__1; 00021 00022 /* Local variables */ 00023 extern logical lsame_(char *, char *); 00024 extern /* Subroutine */ int xerbla_(char *, integer *), dsptrf_( 00025 char *, integer *, doublereal *, integer *, integer *), 00026 dsptrs_(char *, integer *, integer *, doublereal *, integer *, 00027 doublereal *, integer *, integer *); 00028 00029 00030 /* -- LAPACK driver routine (version 3.2) -- */ 00031 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00032 /* November 2006 */ 00033 00034 /* .. Scalar Arguments .. */ 00035 /* .. */ 00036 /* .. Array Arguments .. */ 00037 /* .. */ 00038 00039 /* Purpose */ 00040 /* ======= */ 00041 00042 /* DSPSV computes the solution to a real system of linear equations */ 00043 /* A * X = B, */ 00044 /* where A is an N-by-N symmetric matrix stored in packed format and X */ 00045 /* and B are N-by-NRHS matrices. */ 00046 00047 /* The diagonal pivoting method is used to factor A as */ 00048 /* A = U * D * U**T, if UPLO = 'U', or */ 00049 /* A = L * D * L**T, if UPLO = 'L', */ 00050 /* where U (or L) is a product of permutation and unit upper (lower) */ 00051 /* triangular matrices, D is symmetric and block diagonal with 1-by-1 */ 00052 /* and 2-by-2 diagonal blocks. The factored form of A is then used to */ 00053 /* solve the system of equations A * X = B. */ 00054 00055 /* Arguments */ 00056 /* ========= */ 00057 00058 /* UPLO (input) CHARACTER*1 */ 00059 /* = 'U': Upper triangle of A is stored; */ 00060 /* = 'L': Lower triangle of A is stored. */ 00061 00062 /* N (input) INTEGER */ 00063 /* The number of linear equations, i.e., the order of the */ 00064 /* matrix A. N >= 0. */ 00065 00066 /* NRHS (input) INTEGER */ 00067 /* The number of right hand sides, i.e., the number of columns */ 00068 /* of the matrix B. NRHS >= 0. */ 00069 00070 /* AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2) */ 00071 /* On entry, the upper or lower triangle of the symmetric matrix */ 00072 /* A, packed columnwise in a linear array. The j-th column of A */ 00073 /* is stored in the array AP as follows: */ 00074 /* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; */ 00075 /* if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. */ 00076 /* See below for further details. */ 00077 00078 /* On exit, the block diagonal matrix D and the multipliers used */ 00079 /* to obtain the factor U or L from the factorization */ 00080 /* A = U*D*U**T or A = L*D*L**T as computed by DSPTRF, stored as */ 00081 /* a packed triangular matrix in the same storage format as A. */ 00082 00083 /* IPIV (output) INTEGER array, dimension (N) */ 00084 /* Details of the interchanges and the block structure of D, as */ 00085 /* determined by DSPTRF. If IPIV(k) > 0, then rows and columns */ 00086 /* k and IPIV(k) were interchanged, and D(k,k) is a 1-by-1 */ 00087 /* diagonal block. If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, */ 00088 /* then rows and columns k-1 and -IPIV(k) were interchanged and */ 00089 /* D(k-1:k,k-1:k) is a 2-by-2 diagonal block. If UPLO = 'L' and */ 00090 /* IPIV(k) = IPIV(k+1) < 0, then rows and columns k+1 and */ 00091 /* -IPIV(k) were interchanged and D(k:k+1,k:k+1) is a 2-by-2 */ 00092 /* diagonal block. */ 00093 00094 /* B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) */ 00095 /* On entry, the N-by-NRHS right hand side matrix B. */ 00096 /* On exit, if INFO = 0, the N-by-NRHS solution matrix X. */ 00097 00098 /* LDB (input) INTEGER */ 00099 /* The leading dimension of the array B. LDB >= max(1,N). */ 00100 00101 /* INFO (output) INTEGER */ 00102 /* = 0: successful exit */ 00103 /* < 0: if INFO = -i, the i-th argument had an illegal value */ 00104 /* > 0: if INFO = i, D(i,i) is exactly zero. The factorization */ 00105 /* has been completed, but the block diagonal matrix D is */ 00106 /* exactly singular, so the solution could not be */ 00107 /* computed. */ 00108 00109 /* Further Details */ 00110 /* =============== */ 00111 00112 /* The packed storage scheme is illustrated by the following example */ 00113 /* when N = 4, UPLO = 'U': */ 00114 00115 /* Two-dimensional storage of the symmetric matrix A: */ 00116 00117 /* a11 a12 a13 a14 */ 00118 /* a22 a23 a24 */ 00119 /* a33 a34 (aij = aji) */ 00120 /* a44 */ 00121 00122 /* Packed storage of the upper triangle of A: */ 00123 00124 /* AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ] */ 00125 00126 /* ===================================================================== */ 00127 00128 /* .. External Functions .. */ 00129 /* .. */ 00130 /* .. External Subroutines .. */ 00131 /* .. */ 00132 /* .. Intrinsic Functions .. */ 00133 /* .. */ 00134 /* .. Executable Statements .. */ 00135 00136 /* Test the input parameters. */ 00137 00138 /* Parameter adjustments */ 00139 --ap; 00140 --ipiv; 00141 b_dim1 = *ldb; 00142 b_offset = 1 + b_dim1; 00143 b -= b_offset; 00144 00145 /* Function Body */ 00146 *info = 0; 00147 if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) { 00148 *info = -1; 00149 } else if (*n < 0) { 00150 *info = -2; 00151 } else if (*nrhs < 0) { 00152 *info = -3; 00153 } else if (*ldb < max(1,*n)) { 00154 *info = -7; 00155 } 00156 if (*info != 0) { 00157 i__1 = -(*info); 00158 xerbla_("DSPSV ", &i__1); 00159 return 0; 00160 } 00161 00162 /* Compute the factorization A = U*D*U' or A = L*D*L'. */ 00163 00164 dsptrf_(uplo, n, &ap[1], &ipiv[1], info); 00165 if (*info == 0) { 00166 00167 /* Solve the system A*X = B, overwriting B with X. */ 00168 00169 dsptrs_(uplo, n, nrhs, &ap[1], &ipiv[1], &b[b_offset], ldb, info); 00170 00171 } 00172 return 0; 00173 00174 /* End of DSPSV */ 00175 00176 } /* dspsv_ */