՝ABS C :LINPACK H tSGEDI C 2FBLAS C ;=-SGECO C xMOUT C SGEFA C -$`SGESL C 2HILGEN C 0SLDR C C !CODR C cDIDR C }&D-SGEFA C -SGESL C 2SLDR C C /* Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 */ float Abs ( x ) float x ; { if ( x < 0.0 ) return( -x ) ; else return( x ) ; } /* Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 */ /* Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 */ #define ldx 10 int lda ; int n ; int PrintCode ; float A[ldx][ldx] ; float RCond ; int InfoCode ; int IPvt[ldx] ; float Det[2] ; float b[ldx] ; int JobCode ; float Work[ldx] ; /* Copyright (C) Adam Fritz, 133 Main St., Afton, NY 13730 */ /* Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 */ #define TEN 10.0 #define iTEN 10 SGEDI ( A, lda, n, IPvt, Det, Work, JobCode ) float A[][ldx] ; int lda, n ; int IPvt[] ; float Det[] ; float Work[] ; int *JobCode ; /* Program: */ /* */ /* SGEDI - computes the determinant and inverse of */ /* a matrix using the factors computed by */ /* SGECO or SGEFA. */ /* */ /* Version: Date: */ /* */ /* C 02/22/85 */ /* FORTRAN 08/14/78 */ /* */ /* Description: */ /* */ /* On entry; */ /* */ /* A - Output from SGECO or SGEFA. */ /* lda - Leading dimension of matrix A. */ /* n - Order of matrix A. */ /* IPvt - Pivot index vector from SGECO or */ /* SGEFA. */ /* Work - Work vector. Contents destroyed. */ /* JobCode - */ /* =11 compute determinant and inverse. */ /* =01 compute inverse only. */ /* =10 compute determinant only. */ /* */ /* On return; */ /* */ /* A - Inverse of original matrix if */ /* requested, otherwise unchanged. */ /* Det - Determinant of original matrix if */ /* requested, otherwise ignored; */ /* determinant = Det[1]*10**Det[2] */ /* with; */ /* 1.0 <= abs(Det[1]) < 10.0 */ /* or */ /* Det[1] = 0.0. */ /* */ /* Error condition; */ /* */ /* Division by zero will occur if the input */ /* factor contains a zero on the diagonal and */ /* the inverse is requested. It will not occur */ /* if the subroutines are called correctly and */ /* if SGECO has set RCond > 0.0 or SGEFA has */ /* set InfoCode = 0. */ /* */ /* Subprograms: */ /* */ /* SAXPY from BLAS */ /* SSCAL from BLAS */ /* SSWAP from BLAS */ /* */ /* Authors: */ /* */ /* Cleve Moler */ /* University of New Mexico */ /* Argonne National Laboratories */ /* */ /* C and Pascal translations; */ /* */ /* Adam Fritz */ /* 133 Main Street */ /* Afton, New York 13730 */ { float Abs() ; float d ; float e ; float t ; int i, j, k, l ; int kb, kp1, nm1 ; /* Compute Determinant */ if ( *JobCode/iTEN != 0 ) { d = 1.0 ; e = 0.0 ; i = 0 ; do { if ( IPvt[i] != i ) d = -d ; d = A[i][i] * d ; if ( d != 0.0 ) { while ( Abs(d) < 1.0 ) { d *= TEN ; e -= 1.0 ; } while ( Abs(d) >= 10.0 ) { d /= TEN ; e += 1.0 ; } } } while ((++i < n) && (d != 0.0)) ; Det[0] = d ; Det[1] = e ; } /* Compute inverse(U) */ if ( (*JobCode % iTEN) != 0 ) { for ( k = 0; k < n; k++ ) { A[k][k] = 1.0/A[k][k] ; t = -A[k][k] ; SSCAL(k, t, &A[0][k], lda) ; if ( k < n-1 ) { kp1 = k + 1 ; for ( j = kp1; j < n; j++ ) { t = A[k][j] ; A[k][j] = 0.0 ; SAXPY(k+1, t, &A[0][k], lda, &A[0][j], lda) ; } } } /* Form inverse(U) * inverse(L) */ if ( n > 1 ) { nm1 = n - 1 ; for ( kb = 0; kb < nm1; kb++ ) { k = n - kb - 2 ; kp1 = k + 1 ; for ( i = kp1; i < n; i++ ) { Work[i] = A[i][k] ; A[i][k] = 0.0 ; } for ( j = kp1; j < n; j++ ) { t = Work[j] ; SAXPY(n, t, &A[0][j], lda, &A[0][k], lda) ; } l = IPvt[k] ; if ( l != k ) SSWAP(n, &A[0][k], lda, &A[0][l], lda) ; } } } } /* Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 */ /* Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 */ int ISAMAX ( n, x, incx ) int n ; float x[] ; int incx ; /* Find index of element having maximum absolute value.*/ /* */ /* Jack Dongarra, LINPACK, 3/11/78. */ /* Adam Fritz, C, 2/22/85. */ { float Abs() ; float smax ; int i, ix ; int imax ; imax = 0 ; if ( n > 0 ) { if ( n > 1 ) { if ( incx > 1 ) { /* incx > 1 */ ix = 0 ; smax = Abs(x[0]) ; ix = incx ; for ( i = 1; i < n; i++ ) { if ( Abs(x[ix]) > smax ) { imax = i ; smax = Abs(x[ix]) ; } ix = ix + incx ; } } /* incx = 1 */ else { smax = Abs(x[0]) ; for ( i = 1; i < n; i++ ) if ( Abs(x[i]) > smax ) { imax = i ; smax = Abs(x[i]) ; } } } } return ( imax ) ; } /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ float SASUM ( n, x, incx ) int n ; float x[] ; int incx ; /* Forms sum of absolute values. */ /* */ /* Jack Dongarra, LINPACK, 3/11/78. */ /* Adam Fritz, C, 2/22/85. */ { float Abs() ; float stemp ; int i, ix ; stemp = 0.0 ; if ( n > 0 ) { if ( incx > 1 ) { /* incx > 1 */ ix = 0 ; for ( i = 0; i < n; i++ ) { stemp = stemp + Abs(x[ix]) ; ix = ix + incx ; } } /* incx = 1 */ else for ( i = 0; i < n; i++ ) stemp = stemp + Abs(x[i]) ; } return( stemp ) ; } /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ SAXPY ( n, sa, x, incx, y, incy ) int n ; float sa ; float x[] ; int incx ; float y[] ; int incy ; /* Compute constant times a vector plus a vector. */ /* */ /* Jack Dongarra, LINPACK, 3/11/78. */ /* Adam Fritz, C, 2/22/85. */ { int i, ix, iy ; if ( n > 0 ) { if ( sa != 0.0 ) { if ( (incx != 1) || (incy != 1) ) { /* incx != incy or incx != 1 */ ix = 0 ; iy = 0 ; if ( incx < 0 ) ix = (-n + 1) * incx ; if ( incy < 0 ) iy = (-n + 1) * incy ; for ( i = 0; i < n; i++ ) { y[iy] = y[iy] + sa * x[ix] ; ix = ix + incx ; iy = iy + incy ; } } /* incx, incy = 1 */ else for ( i = 0; i < n; i++ ) y[i] = y[i] + sa * x[i] ; } } } /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ SCOPY ( n, x, incx, y, incy ) int n ; float x[] ; int incx ; float y[] ; int incy ; /* Copies a vector, x, to a vector, y. */ /* */ /* Jack Dongarra, LINPACK, 3/11/78. */ /* Adam Fritz, C, 2/22/85. */ { int i, ix, iy ; if ( n > 0 ) { if ( (incx != 1) || (incy != 1) ) { /* incx, incy != 1 */ ix = 0 ; iy = 0 ; if ( incx < 0 ) ix = (-n + 1) * incx ; if ( incy < 0 ) iy = (-n + 1) * incy ; for ( i = 0; i < n; i++ ) { y[iy] = x[ix] ; ix = ix + incx ; iy = iy + incy ; } } /* incx, incx = 1 */ else for ( i = 0; i < n; i++ ) y[i] = x[i] ; } } /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ float SDOT ( n, x, incx, y, incy ) int n ; float x[] ; int incx ; float y[] ; int incy ; /* Computes dot product of two vectors. */ /* */ /* Jack Dongarra, LINPACK, 3/11/78. */ /* Adam Fritz, C, 2/22/85. */ { float stemp ; int i, ix, iy ; stemp = 0.0 ; if ( n > 0 ) { if ( (incx != 1) || (incy != 1) ) { /* incx != incy or incx != 1 */ ix = 0 ; iy = 0 ; if ( incx < 0 ) ix = (-n + 1) * incx ; if ( incy < 0 ) iy = (-n + 1) * incy ; for ( i = 0; i < n; i++ ) { stemp = stemp + x[ix] * y[iy] ; ix = ix + incx ; iy = iy + incy ; } } /* incx, incy = 1 */ else for ( i = 0; i < n; i++ ) stemp = stemp + x[i] * y[i] ; } return( stemp ) ; } /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ SSCAL ( n, sa, x, incx ) int n ; float sa ; float x[] ; int incx ; /* Computes vector scaled by a constant. */ /* */ /* Jack Dongarra, LINPACK, 3/11/78. */ /* Adam Fritz, C, 2/22/85. */ { int i, ix ; if ( n > 0 ) { if ( incx != 1 ) { /* incx != 1 */ ix = 0 ; for ( i = 0; i < n; i++ ) { x[ix] = sa * x[ix] ; ix = ix + incx ; } } /* incx = 1 */ else for ( i = 0; i < n; i++ ) x[i] = sa * x[i] ; } } /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ SSWAP ( n, x, incx, y, incy ) int n ; float x[] ; int incx ; float y[] ; int incy ; /* Interchanges two vectors. */ /* */ /* Jack Dongarra, LINPACK, 3/11/78. */ /* Adam Fritz, C, 2/22/85. */ { float stemp ; int i, ix, iy ; if ( n > 0 ) { if ( (incx != 1) || (incy != 1) ) { /* incx != incy or incx != 1 */ ix = 0 ; iy = 0 ; if ( incx < 0 ) ix = (-n + 1) * incx ; if ( incy < 0 ) iy = (-n + 1) * incy ; for ( i = 0; i < n; i++ ) { stemp = x[ix] ; x[ix] = y[iy] ; y[iy] = stemp ; ix = ix + incx ; iy = iy + incy ; } } /* incx, incy = 1 */ else for ( i = 0; i < n; i++ ) { stemp = x[i] ; x[i] = y[i] ; y[i] = stemp ; } } } /* Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 */ /* Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 */ SGECO ( A, lda, n, IPvt, RCond, z ) float A[][ldx] ; int lda, n ; int IPvt[] ; float *RCond ; float z[] ; /* Program: */ /* */ /* SGECO - Factors a real matrix by Gaussian */ /* elimination and computes the condition */ /* RCond of the matrix. */ /* */ /* Version: Date: */ /* */ /* C 02/22/85 */ /* FORTRAN 08/14/78 */ /* */ /* Description: */ /* */ /* 1 */ /* RCond = ---------------------------------------- */ /* norm(A) * (Estimate of norm(inverse(A))) */ /* */ /* where; */ /* */ /* Estimate = norm(z)/norm(y) */ /* A * z = y */ /* trans(A) * y = E */ /* */ /* and trans(A) is transpose of A. Components of E */ /* are chosen to cause maximum local growth in the */ /* elements of w where; */ /* */ /* trans(U) * w = E */ /* */ /* Vectors are rescaled to avoid overflow. */ /* */ /* If RCond is not needed, SGEFA is slightly faster.*/ /* */ /* To compute: */ /* */ /* A * x = b , follow SGECO by SGESL. */ /* inverse(A) * c , follow SGECO by SGESL. */ /* determinant(A) , follow SGECO by SGEDI. */ /* inverse(A) , follow SGECO by SGEDI. */ /* */ /* On entry; */ /* */ /* A - Matrix to be factored. */ /* lda - Leading dimension of matrix A. */ /* n - Order of matrix A. */ /* */ /* On exit; */ /* */ /* A - Upper triangular matrix and multi- */ /* pliers used to obtain it. Factoriza- */ /* tion can be written; */ /* A = L * U */ /* where L is a product of permutation */ /* and unit lower triangular matrices */ /* and U is upper triangular. */ /* */ /* IPvt - Pivot index vector. */ /* */ /* RCond - Estimate of reciprocal condition of A.*/ /* For system A * x = B , relative per- */ /* turbations in A and b of size */ /* Epsilon may cause relative perturba- */ /* tions in x of size Epsilon/RCond . */ /* If RCond is so small that the logi- */ /* cal expression; */ /* (1.0 + RCond) = 1.0 */ /* is true, then A may be singular to */ /* workingprecision. In particular, */ /* RCond is zero if exact singularity */ /* is detected or the estimate under- */ /* flows. */ /* */ /* z - Work vector whose contents are */ /* usually unimportant. If A is close to */ /* a singular matrix, then z is an */ /* approximate null vector in the sense */ /* that; */ /* norm(A*z) = RCond*norm(A)*norm(z) */ /* */ /* Subprograms: */ /* */ /* SGEFA from LINPACK */ /* SAXPY from BLAS */ /* SDOT from BLAS */ /* SSCAL from BLAS */ /* SASUM from BLAS */ /* */ /* Authors: */ /* */ /* Cleve Moler */ /* University of New Mexico */ /* Argonne National Laboratories */ /* */ /* C and Pascal Translations; */ /* */ /* Adam Fritz */ /* 133 Main Street */ /* Afton, New York 13730 */ { float Abs() ; float SASUM(), SDOT() ; float ek, t, wk, wkm ; float anorm, s, sm, ynorm ; int j, k, kb, kp1, l ; /* Compute 1-norm of A */ anorm = 0.0 ; for ( j = 0; j < n; j++ ) { s = SASUM(n, &A[0][j], lda) ; if ( anorm < s ) anorm = s ; } /* Factor */ SGEFA(A, lda, n, IPvt, &InfoCode) ; /* Solve trans(U) * w = E */ ek = 1.0 ; for ( j = 0; j < n; j++ ) z[j] = 0.0 ; for ( k = 0; k < n; k++ ) { if ( z[k] > 0.0 ) ek = -Abs(ek) ; else ek = Abs(ek) ; if ( Abs(ek-z[k]) > Abs(A[k][k]) ) { s = Abs(A[k][k])/Abs(ek-z[k]) ; SSCAL(n, s, &z[0], 1) ; ek = s * ek ; } wk = ek - z[k] ; wkm = -ek - z[k] ; s = Abs(wk) ; sm = Abs(wkm) ; if ( A[k][k] != 0.0 ) { wk = wk/A[k][k] ; wkm = wkm/A[k][k] ; } else { wk = 1.0 ; wkm = 1.0 ; } kp1 = k + 1 ; if ( kp1 < n ) { for ( j = kp1; j < n; j++ ) { sm = sm + Abs(z[j] + wkm * A[k][j]) ; z[j] = z[j] + wk * A[k][j] ; s = s + Abs(z[j]) ; } if ( s < sm ) { t = wkm - wk ; wk = wkm ; for ( j = kp1; j < n; j++ ) z[j] = z[j] + t * A[k][j] ; } } z[k] = wk ; } s = 1.0/SASUM(n, &z[0], 1) ; SSCAL(n, s, &z[0], 1) ; /* Solve trans(L) * y = w */ for ( kb = 0; kb < n; kb++ ) { k = n - kb - 1 ; if ( k < n-1 ) z[k] = z[k] + SDOT(n-k-1, &A[k+1][k], lda, &z[k+1], 1) ; if ( Abs(z[k]) > 1.0 ) { s = 1.0/Abs(z[k]) ; SSCAL(n, s, &z[0], 1) ; } l = IPvt[k] ; t = z[l] ; z[l] = z[k] ; z[k] = t ; } s = 1.0/SASUM(n, &z[0], 1) ; SSCAL(n, s, &z[0], 1) ; /* Solve L * v = y */ ynorm = 1.0 ; for ( k = 0; k < n; k++ ) { l = IPvt[k] ; t = z[l] ; z[l] = z[k] ; z[k] = t ; if ( k < n-1 ) SAXPY(n-k-1, t, &A[k+1][k], lda, &z[k+1], 1) ; if ( Abs(z[k]) > 1.0 ) { s = 1.0/Abs(z[k]) ; SSCAL(n, s, &z[0], 1) ; ynorm = s * ynorm ; } } s = 1.0/SASUM(n, &z[0], 1) ; SSCAL(n, s, &z[0], 1) ; ynorm = s * ynorm ; /* Solve U * z = v */ for ( kb = 0; kb < n; kb++ ) { k = n - kb - 1 ; if ( Abs(z[k]) > Abs(A[k][k]) ) { s = Abs(A[k][k])/Abs(z[k]) ; SSCAL(n, s, &z[0], 1) ; ynorm = s * ynorm ; } if ( A[k][k] != 0.0 ) z[k] = z[k]/A[k][k] ; else z[k] = 1.0 ; t = -z[k] ; SAXPY(k, t, &A[0][k], lda, &z[0], 1) ; } /* Set znorm = 1.0 */ s = 1.0/SASUM(n, &z[0], 1) ; SSCAL(n, s, &z[0], 1) ; ynorm = s * ynorm ; if ( anorm != 0.0 ) *RCond = ynorm/anorm ; else *RCond = 0.0 ; } /* Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 */ /* Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 */ OUT ( A, lda, n , m ) float A[] ; int lda, n, m ; /* General purpose output routine for the matrix A */ /* which has leading dimension lda and is n by m. */ /* */ /* J.J. Dongarra, SICE, ... (?) */ /* Adam Fritz, C, 2/22/85 */ { int i, j, k ; int ic, icb, ice ; if ( n > 0 ) { if ( m > 1 ) { ic = (m+4)/5 ; icb = 0 ; ice = 5 ; for ( k = 0; k < ic; k++ ) { if (ice > m ) ice = m ; for ( i = 0; i < n; i++ ) { for ( j = icb; j < ice; j++ ) printf("%7e ", A[i*lda+j]) ; printf("\n") ; } ; icb = icb + 5 ; ice = ice + 5 ; printf("\n") ; } } else { ic = (n+4)/5 ; icb = 0 ; ice = 5 ; for ( k = 0; k < ic; k++ ) { if (ice > n ) ice = n ; for ( i = icb; i < ice; i++ ) printf("%7e ", A[i]) ; printf("\n") ; icb = icb + 5 ; ice = ice + 5 ; } ; printf("\n") ; } } } /* Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 */ /* Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 */ SGEFA ( A, lda, n, IPvt, InfoCode ) float A[][ldx] ; int lda, n ; int IPvt[] ; int *InfoCode ; /* Program: */ /* */ /* SGEFA - Factors Real Single Precision Matrix */ /* Using Gaussian Elimination. */ /* */ /* Version: DATE: */ /* */ /* C 02/22/85 */ /* FORTRAN 08/14/78 */ /* */ /* Description: */ /* */ /* SGEFA is usually called by SGECO, but can be */ /* called directly with a saving in time if RCOND */ /* is not needed. */ /* */ /* (time for SGECO) = (1 + 9/n)*(time for SGEFA) */ /* */ /* On entry; */ /* */ /* A - Matrix to be factored. */ /* lda - Leading dimension of matrix A. */ /* n - Order of matrix A. */ /* */ /* On exit; */ /* */ /* A - Upper triangular matrix and multi- */ /* pliers used to obtain it. Factoriza- */ /* tion can be written; */ /* A = L * U */ /* where L is a product of permutation */ /* and unit lower triangular matrices */ /* and U is upper triangular. */ /* */ /* IPvt - Pivot index vector. */ /* */ /* InfoCode - */ /* = 0 Normal value. */ /* = k If u[k,k] = 0.0. This is not an */ /* error condition for this procedure, */ /* but indicates that SGESL or SGEDI */ /* will divide by zero if called. Use */ /* RCOND from SGECO for reliable indi- */ /* cation of singularity. */ /* */ /* Subprograms: */ /* */ /* SAXPY from BLAS */ /* SSCAL from BLAS */ /* ISAMAX from BLAS */ /* */ /* Authors: */ /* */ /* Cleve Moler */ /* University of New Mexico */ /* Argonne National Laboratories */ /* */ /* C and Pascal translations; */ /* */ /* Adam Fritz */ /* 133 Main Street */ /* Afton, New York 13730 */ { int ISAMAX() ; int j, k, l ; int kp1, nm1 ; float t ; /* Gaussian Elimination with Partial Pivoting */ *InfoCode = 0 ; nm1 = n - 1 ; if ( nm1 > 0 ) { for ( k = 0; k < nm1; k++ ) { kp1 = k + 1 ; /* Find l = Pivot Index */ l = ISAMAX(n-k, &A[k][k], lda) + k ; IPvt[k] = l ; /* Zero Pivot Implies Column Triangularized */ if ( A[l][k] != 0.0 ) { /* Interchange If Necessary */ if ( l != k ) { t = A[l][k] ; A[l][k] = A[k][k] ; A[k][k] = t ; } /* Compute Multipliers */ t = -1.0/A[k][k] ; SSCAL(n-k-1, t, &A[k+1][k], lda) ; /* Row Elimination with Column Indexing */ for ( j = kp1; j < n; j++ ) { t = A[l][j] ; if ( l != k ) { A[l][j] = A[k][j] ; A[k][j] = t ; } SAXPY(n-k-1, t, &A[k+1][k], lda, &A[k+1][j], lda) ; } } else *InfoCode = k ; } } IPvt[nm1] = nm1 ; if ( A[nm1][nm1] == 0.0 ) *InfoCode = nm1 ; } /* Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 */ /* Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 */ SGESL ( A, lda, n, IPvt, b, JobCode ) float A[][ldx] ; int lda, n ; int IPvt[] ; float b[] ; int *JobCode ; /* Program: */ /* */ /* SGESL - Solve General Real Linear System. */ /* */ /* Version: Date: */ /* */ /* C 02/22/85 */ /* FORTRAN 08/14/78 */ /* */ /* Description: */ /* */ /* Solves real single precision system; */ /* */ /* A * x = b or trans(A) * x = b */ /* */ /* using factors computed by SGECO or SGEFA. */ /* */ /* On entry; */ /* */ /* A - Output from SGECO or SGEFA */ /* lda - Leading dimension of matrix A */ /* n - Order of matrix A. */ /* IPvt - Pivot index vector from SGECO or */ /* SGEFA. */ /* b - Right hand side vector */ /* JobCode - */ /* =0 solve A * x = b , */ /* >0 solve trans(A) * x = b where */ /* trans(A) is the transpose. */ /* */ /* On exit; */ /* */ /* b - Solution vector x . */ /* */ /* Error Condition; */ /* */ /* Division by zero will occur if the input */ /* factor contains a zero on the diagonal. This */ /* indicates singularity but can be caused by */ /* improper arguments or improper setting of lda.*/ /* It will occur if the procedures are called */ /* correctly and SGECO has set RCOND greater */ /* than 0.0 or SGEFA has set InfoCode */ /* different from 0. */ /* */ /* To compute; */ /* */ /* inverse(A) * c */ /* */ /* where c is a matrix with p columns; */ /* */ /* SGECO(A, lda, n, IPvt, RCOND, z) ; */ /* if RCOND ... ok */ /* for j = 1 to p do */ /* SGESL(A, lda, n, IPvt, c[1,j], false) ; */ /* */ /* Subprograms: */ /* */ /* SAXPY from BLAS */ /* SDOT from BLAS */ /* */ /* Authors: */ /* */ /* Cleve Moler */ /* University of New Mexico */ /* Argonne National Laboratories */ /* */ /* C and Pascal translations; */ /* */ /* Adam Fritz */ /* 133 Main Street */ /* Afton, New York 13730 */ { float SDOT() ; int k, kb, l, nm1 ; float t ; nm1 = n - 1 ; if ( *JobCode == 0 ) { /* Solve A * x = b */ if ( nm1 > 0 ) { for ( k = 0; k < nm1; k++ ) { l = IPvt[k] ; t = b[l] ; if ( l != k ) { b[l] = b[k] ; b[k] = t ; } SAXPY(n-k-1, t, &A[k+1][k], lda, &b[k+1], 1) ; } } /* Solve U * x = y */ for ( kb = 0; kb < n; kb++ ) { k = n - kb - 1 ; b[k] = b[k]/A[k][k] ; t = -b[k] ; SAXPY(k, t, &A[0][k], lda, &b[0], 1) ; } } /* Solve trans(A) * x = b */ else { /* Solve trans(U) * y = b */ for ( k = 0; k < n; k++ ) { t = SDOT(k, &A[0][k], lda, &b[0], 1) ; b[k] = (b[k] - t)/A[k][k] ; } /* Solve trans(L) * x = y */ if ( nm1 > 0 ) { for ( kb = 0; kb < nm1; kb++ ) { k = n - kb - 2 ; b[k] = b[k] + SDOT(n-k-1, &A[k+1][k], lda, &b[k+1], 1) ; l = IPvt[k] ; if ( l != k ) { t = b[l] ; b[l] = b[k] ; b[k] = t ; } } } } } /* Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 */ /* Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 */ SYSGEN ( A, lda, n, b ) float A[][ldx], b[] ; int lda, n ; /* Program: */ /* */ /* SYSGEN */ /* */ /* Version: Date: */ /* */ /* C 02/22/85 */ /* */ /* Description: */ /* */ /* Generate a test matrix and test vector using */ /* Hilbert coefficients and set vector to matrix */ /* diagonal elements. The matrix is symmetric and */ /* poorly conditioned. */ /* */ /* Author: */ /* */ /* Adam Fritz */ /* 133 Main Street */ /* Afton, New York 13730 */ { int i, j, ip1 ; float t ; /* Validate Leading Dimension */ if ( lda > 0 ) { /* Validate Order */ if ((n > 1) && (n <= lda)) { /* Form Matrix and Vector */ for ( i = 0; i < n; i++ ) { t = 1.0/(2*i+1) ; A[i][i] = t ; b[i] = t ; ip1 = i + 1 ; for ( j = ip1; j < n; j++ ) { t = 1.0/(i+j+1) ; A[i][j] = t ; A[j][i] = t ; } } } else { printf("Error: Invalid MATRIX Order, %d\n", n) ; bios(0) ; } } else { printf("Error: Invalid LEADING DIMENSION, %d\n", lda) ; bios(0) ; } } /* Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 */ /* Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 */ #include "stdio.h" #include "linpack.h" #include "Abs.c" #include "Hilgen.c" /* Test System Generator */ #include "BLAS.c" /* Basic Linear Algebra */ #include "SGEFA.c" /* LINPACK Factor */ #include "SGESL.c" /* LINPACK Solve */ #include "OUT.c" /* SICE Output Routine */ main () /* Program: */ /* */ /* LINPACK SGEFA and SGESL Test Driver. */ /* */ /* Version: Date: */ /* */ /* C 02/22/85 */ /* */ /* Description: */ /* */ /* Uses LINPACK SGEFA and SGESL to solve; */ /* */ /* A * x = b */ /* */ /* where; */ /* */ /* A is set up as a Hilbert matrix and b is a */ /* 'Hilbert' vector of specified order using */ /* procedure SYSGEN, x is computed using SGEFA */ /* to factor a, SGESL solves the linear system */ /* in single precision, and results are shown */ /* conditioned on a print code. */ /* */ /* Author: */ /* */ /* Adam Fritz */ /* 133 Main Street */ /* Afton, New York 13730 */ { /* Initialize */ lda = ldx ; /* Announcement */ printf("LINPACK SGEFA and SGESL Test Program.\n") ; printf("\n") ; /* Get Order */ n = 0 ; while ((n < 1) || (n > lda)) { printf("Order: ") ; scanf("%d", &n) ; } /* Get Print Code */ printf("Print Code: ") ; scanf("%d", &PrintCode) ; /* Generate Test System */ SYSGEN(A, lda, n, b) ; if ( PrintCode > 0 ) { printf("\n") ; printf("Original System (by column):\n") ; printf("\n") ; OUT(A, lda, n, n) ; OUT(b, lda, n, 1) ; } /* Factor the Matrix */ printf("begin: SGEFA\n") ; SGEFA(A, lda, n, IPvt, &InfoCode) ; printf("end: SGEFA\n") ; if ( PrintCode > 0 ) { printf("\n") ; printf("Factored Matrix (by column):\n") ; printf("\n") ; OUT(A, lda, n, n) ; } /* Solve the Linear System */ if ( InfoCode == 0 ) { JobCode = 0 ; printf("begin: SGESL\n") ; SGESL(A, lda, n, IPvt, b, &JobCode) ; printf("end: SGESL\n") ; if ( PrintCode > 0 ) { printf("\n") ; printf("Solution Vector:\n") ; printf("\n") ; OUT(b, lda, n, 1) ; } } /* Zero Diagonal Element Produced */ else printf("Error: Possible Singular System.\n") ; /* Done */ printf("End of Test.\n") ; } /* Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 */ /* Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 */ #include "stdio.h" #include "linpack.h" #include "Abs.c" #include "Hilgen.c" /* Test System Generator */ #include "BLAS.c" /* Basic Linear Algebra */ #include "SGEFA.c" /* LINPACK Factor */ #include "SGECO.c" /* LINPACK Condition */ #include "OUT.c" /* SICE Output Routine */ main () /* Program: */ /* */ /* LINPACK SGECO and SGEFA Test Driver. */ /* */ /* Version: Date: */ /* */ /* C 02/22/85 */ /* */ /* Description: */ /* */ /* Uses LINPACK SGECO and SGEFA to compute */ /* condition estimate RCond for matrices. */ /* A is set up as a Hilbert matrix of specified */ /* order and SGECO is called to compute the RCond */ /* measure. If PrintCode != 0 then printout will */ /* include RCond and 'folded' RCond where; */ /* folded RCond = (1.0+RCond)-1.0 */ /* */ /* Author: */ /* */ /* Adam Fritz */ /* 133 Main Street */ /* Afton, New York 13730 */ /* */ { /* Initialize */ lda = ldx ; /* Announcement */ printf("LINPACK SGECO and SGEFA Test Program.\n") ; printf("\n") ; /* Get Order */ n = 0 ; while ((n < 1) || (n > lda)) { printf("Order: ") ; scanf("%d", &n) ; } /* Get Print Code */ printf("Print Code: ") ; scanf("%d", &PrintCode) ; /* Generate Test System */ SYSGEN(A, lda, n, b) ; if ( PrintCode > 0 ) { printf("\n") ; printf("Original System (by column):\n") ; printf("\n") ; OUT(A, lda, n, n) ; OUT(b, lda, n, 1) ; } /* Compute the Condition */ SGECO(A, lda, n, IPvt, &RCond, Work) ; if ( PrintCode > 0 ) { printf("\n") ; printf("RCond: %7e", RCond) ; RCond += 1.0 ; RCond -= 1.0 ; printf(", RCond: %7e\n", RCond) ; printf("\n") ; } /* Done */ printf("End of Test.\n") ; } /* Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 */ /* Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 */ #include "stdio.h" #include "linpack.h" #include "Abs.c" #include "Hilgen.c" /* Test System Generator */ #include "BLAS.c" /* Basic Linear Algebra */ #include "SGEFA.c" /* LINPACK Factor */ #include "SGECO.c" /* LINPACK Condition */ #include "SGEDI.c" /* LINPACK Determinant and */ #include "OUT.c" /* SICE Output Routine */ main () /* Program: */ /* */ /* LINPACK SGECO and SGEDI Test Driver. */ /* */ /* Version: Date: */ /* */ /* C 02/22/85 */ /* */ /* Description: */ /* */ /* Uses LINPACK SGECO to compute RCond for matrix */ /* A and LINPACK SGEDI to compute determinant and */ /* inverse of matrix A where A is set up as a */ /* Hilbert matrix of specified order using SYSGEN. */ /* Uses SGEFA and SGEDI to invert inverse to com- */ /* with original. */ /* */ /* Author: */ /* */ /* Adam Fritz */ /* 133 Main Street */ /* Afton, New York 13730 */ /* */ { float SDOT() ; int i, j ; float AA[ldx][ldx] ; float P[ldx][ldx] ; /* Initialize */ lda = ldx ; /* Announcement */ printf("LINPACK SGECO and SDEDI Test Program.") ; printf("\n") ; /* Get Order */ n = 0 ; while ((n < 1) || (n > lda)) { printf("Order: ") ; scanf("%d", &n) ; } /* Get Print Code */ printf("Print Code: ") ; scanf("%d", &PrintCode) ; /* Generate Test System */ SYSGEN(A, lda, n, b) ; if ( PrintCode > 0 ) { printf("\n") ; printf("Original System (by column):\n") ; printf("\n") ; OUT(A, lda, n, n) ; OUT(b, lda, n, 1) ; } /* Copy Original by Rows */ for ( i = 0; i < n; i++ ) SCOPY(n, &A[i][0], 1, &AA[i][0], 1) ; /* Compute Condition */ SGECO(A, lda, n, IPvt, &RCond, Work) ; if ( PrintCode > 0 ) { printf("\n") ; printf("RCond: %7e",RCond) ; RCond += 1.0 ; RCond -= 1.0 ; printf(", RCond: %7e\n", RCond) ; printf("\n") ; } /* Compute Determinant and Inverse */ if ( RCond != 0.0 ) { JobCode = 11 ; SGEDI(A, lda, n, IPvt, Det, Work, &JobCode) ; if ( PrintCode > 0 ) { printf("\n") ; printf("Inverse (by column)\n:") ; printf("\n") ; OUT(A, lda, n, n) ; printf("Determinant: %7fe%-3d\n", Det[0], (int) Det[1]) ; } /* Original Times Inverse */ for ( i = 0; i < n; i++ ) for ( j = 0; j < n; j++ ) P[i][j] = SDOT(n, &AA[i][0], 1, &A[0][j], lda) ; if ( PrintCode > 0 ) { printf("\n") ; printf("Original times Inverse (by column):\n") ; printf("\n") ; OUT(P, lda, n, n) ; } /* Try to Restore Original */ SGEFA(A, lda, n, IPvt, &InfoCode) ; if ( InfoCode == 0 ) { SGEDI(A, lda, n, IPvt, Det, Work, &JobCode) ; if ( PrintCode > 0 ) { printf("Inverse of Inverse (by column):\n") ; printf("\n") ; OUT(A, lda, n, n) ; } } else printf("Error: Inverse is Singular.\n") ; } /* Zero Diagonal Element Produced */ else printf("Error: Original is Singular.\n") ; /* Done */ printf("End of Test.\n") ; } /* Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 */