LINPACK CON9LINPACK TYP8LINPACK VAR #~HILGEN PASOUT PAS!qSGECO PAS/MBLAS PAS|MUWSGEFA PAS-tSGESL PAS2ASGEDI PAS(4U/SLDR PAS\#CODR PASz[DIDR PAS$G t8ßqpjpÏpqâqtB $% ? (Y/N): {Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730} const lda = 10 ; {Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730} {Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730} type MATRIX = array[1..lda,1..lda] of real ; VECTOR = array[1..lda] of real ; IVECTOR = array[1..lda] of integer ; DVECTOR = array[1..2] of real ; RowPointer = ^Row ; Row = record s : VECTOR end ; {Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730} {Copyright (C) Adam Fritz, 133 Main St., Afton, NY 13730} var n : integer ; PrintCode : integer ; A : MATRIX ; RCond : real ; InfoCode : integer ; IPvt : IVECTOR ; Det : DVECTOR ; b : VECTOR ; JobCode : integer ; Work : VECTOR ; {Copyright (C) Adam Fritz, 133 Main St., Afton, NY 13730} {Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730} procedure SYSGEN ( var a : MATRIX ; lda, n : integer ; var b : VECTOR ) ; { } { PROGRAM: SYSGEN } { } { VERSION: 1.0/TURBO DATE: 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 } { } var i, j, nm1, ip1 : integer ; t : real ; begin { Validate Leading Dimension } if lda > 0 then begin { Validate Order } if (n > 1) and (n <= lda) then begin { Form Matrix and Vector } nm1 := n - 1 ; if n > 1 then begin for i := 1 to nm1 do begin t := 1.0/(2*i-1) ; a[i,i] := t ; b[i] := t ; ip1 := i + 1 ; for j := ip1 to n do begin t := 1.0/(i+j-1) ; a[i,j] := t ; a[j,i] := t end end end ; t := 1.0/(2*n-1) ; a[n,n] := t ; b[n] := t end else begin writeln('Error: Invalid MATRIX Order, ', n,'.') ; BIOS(0) end end else begin writeln('Error: Invalid LEADING DIMENSION, ', lda) ; BIOS(0) end end ; {Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730} {Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730} procedure OUT ( var sa : real ; lda, n , m : integer ) ; { } { General purpose output routine for the matrix A } { which has leading dimension lda and is n by m. } { } { J.J. Dongarra, SICE, ... (?) } { Adam Fritz, TURBO Pascal, 2/22/85 } { } var i, j, k : integer ; ic, icb, ice : integer ; a : RowPointer ; begin if n > 0 then begin a := Ptr(Addr(sa)) ; if m > 1 then begin ic := (m + 4) div 5 ; icb := 1 ; ice := 5 ; for k := 1 to ic do begin if ice > m then ice := m ; for i := 1 to n do begin for j := icb to ice do write (a^.s[(i-1)*lda+j]:14, ' ') ; writeln end ; icb := icb + 5 ; ice := ice + 5 ; writeln end end else begin ic := (n + 4) div 5 ; icb := 1 ; ice := 5 ; for k := 1 to ic do begin if ice > n then ice := n ; for i := icb to ice do write (a^.s[i]:14, ' ') ; writeln ; icb := icb + 5 ; ice := ice + 5 end ; writeln end end end ; {Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730} {Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730} procedure SGECO ( var A : MATRIX ; lda, n : integer ; var IPvt : IVECTOR ; var RCond : real ; var z : VECTOR ) ; { } { PROGRAM: } { } { SGECO - Factors a real matrix by Gaussian } { elimination and computes the condition } { RCond of the matrix. } { } { VERSION: DATE: } { } { 1.0/TURBO Pascal 2.0 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 } { } { PASCAL Translation; } { } { Adam Fritz } { 133 Main Street } { Afton, New York 13730 } { } var ek, t, wk, wkm : real ; anorm, s, sm, ynorm : real ; InfoCode : integer ; j, k, kb, kp1, l : integer ; begin { Compute 1-norm of A } anorm := 0.0 ; for j := 1 to n do begin s := SASUM(n, A[1,j], lda) ; if anorm < s then anorm := s end ; { Factor } SGEFA(A, lda, n, IPvt, InfoCode) ; { Solve trans(U) * w = E } ek := 1.0 ; for j := 1 to n do z[j] := 0.0 ; for k := 1 to n do begin if z[k] > 0.0 then ek := -Abs(ek) else ek := Abs(ek) ; if Abs(ek-z[k]) > Abs(A[k,k]) then begin s := Abs(A[k,k])/Abs(ek-z[k]) ; SSCAL(n, s, z[1], 1) ; ek := s * ek end ; wk := ek - z[k] ; wkm := -ek - z[k] ; s := Abs(wk) ; sm := Abs(wkm) ; if A[k,k] <> 0.0 then begin wk := wk/A[k,k] ; wkm := wkm/A[k,k] end else begin wk := 1.0 ; wkm := 1.0 end ; kp1 := k + 1 ; if kp1 <= n then begin for j := kp1 to n do begin sm := sm + Abs(z[j] + wkm * A[k,j]) ; z[j] := z[j] + wk * A[k,j] ; s := s + Abs(z[j]) end ; if s < sm then begin t := wkm - wk ; wk := wkm ; for j := kp1 to n do z[j] := z[j] + t * A[k,j] end end ; z[k] := wk end ; s := 1.0/SASUM(n, z[1], 1) ; SSCAL(n, s, z[1], 1) ; { Solve trans(L) * y = w } for kb := 1 to n do begin k := n + 1 - kb ; if k < n then z[k] := z[k] + SDOT(n-k, A[k+1,k], lda, z[k+1], 1) ; if Abs(z[k]) > 1.0 then begin s := 1.0/Abs(z[k]) ; SSCAL(n, s, z[1], 1) end ; l := IPvt[k] ; t := z[l] ; z[l] := z[k] ; z[k] := t end ; s := 1.0/SASUM(n, z[1], 1) ; SSCAL(n, s, z[1], 1) ; { Solve L * v = y } ynorm := 1.0 ; for k := 1 to n do begin l := IPvt[k] ; t := z[l] ; z[l] := z[k] ; z[k] := t ; if k < n then SAXPY(n-k, t, A[k+1,k], lda, z[k+1], 1) ; if Abs(z[k]) > 1.0 then begin s := 1.0/Abs(z[k]) ; SSCAL(n, s, z[1], 1) ; ynorm := s * ynorm end end ; s := 1.0/SASUM(n, z[1], 1) ; SSCAL(n, s, z[1], 1) ; ynorm := s * ynorm ; { Solve U * z = v } for kb := 1 to n do begin k := n + 1 - kb ; if Abs(z[k]) > Abs(A[k,k]) then begin s := Abs(A[k,k])/Abs(z[k]) ; SSCAL(n, s, z[1], 1) ; ynorm := s * ynorm end ; if A[k,k] <> 0.0 then z[k] := z[k]/A[k,k] else z[k] := 1.0 ; t := -z[k] ; SAXPY(k-1, t, A[1,k], lda, z[1], 1) end ; { Set znorm = 1.0 } s := 1.0/SASUM(n, z[1], 1) ; SSCAL(n, s, z[1], 1) ; ynorm := s * ynorm ; if anorm <> 0.0 then RCond := ynorm/anorm else RCond := 0.0 end ; {Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730} {Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 } function ISAMAX ( n : integer ; var sx : real ; incx : integer ) : integer ; { } { Find index of element having maximum absolute value.} { } { Jack Dongarra, LINPACK, 3/11/78. } { Adam Fritz, TURBO Pascal, 2/22/85. } { } var smax : real ; i, ix : integer ; x : RowPointer ; begin ISAMAX := 0 ; if n > 0 then begin ISAMAX := 1 ; if n > 1 then begin x := Ptr(Addr(sx)) ; if incx > 1 then begin { incx > 1 } ix := 1 ; smax := Abs(x^.s[1]) ; ix := ix + incx ; for i := 2 to n do begin if Abs(x^.s[ix]) > smax then begin ISAMAX := i ; smax := Abs(x^.s[ix]) end ; ix := ix + incx end end { incx = 1 } else begin smax := Abs(x^.s[1]) ; for i := 2 to n do if Abs(x^.s[i]) > smax then begin ISAMAX := i ; smax := Abs(x^.s[i]) end end end end end ; {~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~} function SASUM ( n : integer ; var sx : real ; incx : integer ) : real ; { } { Forms sum of absolute values. } { } { Jack Dongarra, LINPACK, 3/11/78. } { Adam Fritz, TURBO Pascal, 2/22/85. } { } var stemp : real ; i, ix : integer ; x : RowPointer ; begin stemp := 0.0 ; if n > 0 then begin x := Ptr(Addr(sx)) ; if incx > 1 then begin { incx > 1 } ix := 1 ; for i := 1 to n do begin stemp := stemp + Abs(x^.s[ix]) ; ix := ix + incx end end { incx = 1 } else for i := 1 to n do stemp := stemp + Abs(x^.s[i]) end ; SASUM := stemp end ; {~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~} procedure SAXPY ( n : integer ; sa : real ; var sx : real ; incx : integer ; var sy : real ; incy : integer ) ; { } { Compute constant times a vector plus a vector. } { } { Jack Dongarra, LINPACK, 3/11/78. } { Adam Fritz, TURBO Pascal, 2/22/85. } { } var i, ix, iy : integer ; x, y : RowPointer ; begin if n > 0 then begin if sa <> 0.0 then begin x := Ptr(Addr(sx)) ; y := Ptr(Addr(sy)) ; if (incx <> 1) or (incy <> 1) then begin { incx <> incy or incx <> 1 } ix := 1 ; iy := 1 ; if incx < 0 then ix := (-n + 1) * incx + 1 ; if incy < 0 then iy := (-n + 1) * incy + 1 ; for i := 1 to n do begin y^.s[iy] := y^.s[iy] + sa * x^.s[ix] ; ix := ix + incx ; iy := iy + incy end end { incx, incy = 1 } else for i := 1 to n do y^.s[i] := y^.s[i] + sa * x^.s[i] end end end ; {~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~} procedure SCOPY ( n : integer ; var sx : real ; incx : integer ; var sy : real ; incy : integer ) ; { } { Copies a vector, x, to a vector, y. } { } { Jack Dongarra, LINPACK, 3/11/78. } { Adam Fritz, TURBO Pascal, 2/22/85. } { } var i, ix, iy : integer ; x, y : RowPointer ; begin if n > 0 then begin x := Ptr(Addr(sx)) ; y := Ptr(Addr(sy)) ; if (incx <> 1) or (incy <> 1) then begin { incx, incy <> 1 } ix := 1 ; iy := 1 ; if incx < 0 then ix := (-n + 1) * incx + 1 ; if incy < 0 then iy := (-n + 1) * incy + 1 ; for i := 1 to n do begin y^.s[iy] := x^.s[ix] ; ix := ix + incx ; iy := iy + incy end end { incx, incx = 1 } else for i := 1 to n do y^.s[i] := x^.s[i] end end ; {~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~} function SDOT ( n : integer ; var sx : real ; incx : integer ; var sy : real ; incy : integer ) : real ; { } { Computes dot product of two vectors. } { } { Jack Dongarra, LINPACK, 3/11/78. } { Adam Fritz, TURBO Pascal, 2/22/85. } { } var stemp : real ; i, ix, iy : integer ; x, y : RowPointer ; begin stemp := 0.0 ; if n > 0 then begin x := Ptr(Addr(sx)) ; y := Ptr(Addr(sy)) ; if (incx <> 1) or (incy <> 1) then begin { incx <> incy or incx <> 1 } ix := 1 ; iy := 1 ; if incx < 0 then ix := (-n + 1) * incx + 1 ; if incy < 0 then iy := (-n + 1) * incy + 1 ; for i := 1 to n do begin stemp := stemp + x^.s[ix] * y^.s[iy] ; ix := ix + incx ; iy := iy + incy end end { incx, incy = 1 } else for i := 1 to n do stemp := stemp + x^.s[i] * y^.s[i] end ; SDOT := stemp end ; {~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~} procedure SSCAL ( n : integer ; sa : real ; var sx : real ; incx : integer ) ; { } { Computes vector scaled by a constant. } { } { Jack Dongarra, LINPACK, 3/11/78. } { Adam Fritz, TURBO Pascal, 2/22/85. } { } var i, ix : integer ; x : RowPointer ; begin if n > 0 then begin x := Ptr(Addr(sx)) ; if incx <> 1 then begin { incx <> 1 } ix := 1 ; for i := 1 to n do begin x^.s[ix] := sa * x^.s[ix] ; ix := ix + incx end end { incx = 1 } else for i := 1 to n do x^.s[i] := sa * x^.s[i] end end ; {~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~} procedure SSWAP ( n : integer ; var sx : real ; incx : integer ; var sy : real ; incy : integer ) ; { } { Interchanges two vectors. } { } { Jack Dongarra, LINPACK, 3/11/78. } { Adam Fritz, TURBO Pascal, 2/22/85. } { } var stemp : real ; i, ix, iy : integer ; x, y : RowPointer ; begin if n > 0 then begin x := Ptr(Addr(sx)) ; y := Ptr(Addr(sy)) ; if (incx <> 1) or (incy <> 1) then begin { incx <> incy or incx <> 1 } ix := 1 ; iy := 1 ; if incx < 0 then ix := (-n + 1) * incx + 1 ; if incy < 0 then iy := (-n + 1) * incy + 1 ; for i := 1 to n do begin stemp := x^.s[ix] ; x^.s[ix] := y^.s[iy] ; y^.s[iy] := stemp ; ix := ix + incx ; iy := iy + incy end end { incx, incy = 1 } else for i := 1 to n do begin stemp := x^.s[i] ; x^.s[i] := y^.s[i] ; y^.s[i] := stemp end end end ; {Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 } {Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730} procedure SGEFA ( var A : MATRIX ; lda, n : integer ; var IPvt : IVECTOR ; var InfoCode : integer ) ; { } { PROGRAM: } { } { SGEFA - Factors Real Single Precision Matrix } { Using Gaussian Elimination. } { } { VERSION: DATE: } { } { 1.0/TURBO Pascal 2.0 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 } { } { PASCAL translation; } { } { Adam Fritz } { 133 Main Street } { Afton, New York 13730 } { } var j, k, l : integer ; kp1, nm1 : integer ; t : real ; begin { Gaussian Elimination with Partial Pivoting } InfoCode := 0 ; if n > 0 then begin nm1 := n - 1 ; for k := 1 to nm1 do begin kp1 := k + 1 ; { Find l = Pivot Index } l := ISAMAX(n-k+1, A[k,k], lda) + k - 1 ; IPvt[k] := l ; { Zero Pivot Implies Column Triangularized } if A[l,k] <> 0.0 then begin { Interchange if Necessary } if l <> k then begin t := A[l,k] ; A[l,k] := A[k,k] ; A[k,k] := t end ; { Compute Multipliers } t := -1.0/A[k,k] ; SSCAL(n-k, t, A[k+1,k], lda) ; { Row Elimination with Column Indexing } for j := kp1 to n do begin t := A[l,j] ; if l <> k then begin A[l,j] := A[k,j] ; A[k,j] := t end ; SAXPY(n-k, t, A[k+1,k], lda, A[k+1,j], lda) end end else InfoCode := k end end ; IPvt[n] := n ; if A[n,n] = 0.0 then InfoCode := n end ; {Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730} {Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730} procedure SGESL ( var A : MATRIX ; lda, n : integer ; var IPvt : IVECTOR ; var b : VECTOR ; JobCode : integer ) ; { } { PROGRAM: } { } { SGESL - Solve General Real Linear System. } { } { VERSION: DATE: } { } { 1.0/TURBO Pascal 2.0 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 then } { 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 } { } { PASCAL translation; } { } { Adam Fritz } { 133 Main Street } { Afton, New York 13730 } { } var k, kb, l, nm1 : integer ; t : real ; begin nm1 := n - 1 ; if JobCode = 0 then begin { Solve A * x = b } if nm1 > 0 then begin for k := 1 to nm1 do begin l := IPvt[k] ; t := b[l] ; if l <> k then begin b[l] := b[k] ; b[k] := t end ; SAXPY(n-k, t, A[k+1,k], lda, b[k+1], 1) end end ; { Solve U * x = y } for kb := 1 to n do begin k := n + 1 - kb ; b[k] := b[k]/A[k,k] ; t := -b[k] ; SAXPY(k-1, t, A[1,k], lda, b[1], 1) end end { Solve trans(A) * x = b } else begin { Solve trans(U) * y = b } for k := 1 to n do begin t := SDOT(k-1, A[1,k], lda, b[1], 1) ; b[k] := (b[k] - t)/A[k,k] end ; { Solve trans(L) * x = y } if nm1 > 0 then begin for kb := 1 to nm1 do begin k := n - kb ; b[k] := b[k] + SDOT(n-k, A[k+1,k], lda, b[k+1], 1) ; l := IPvt[k] ; if l <> k then begin t := b[l] ; b[l] := b[k] ; b[k] := t end end end end end ; {Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730} {Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730} procedure SGEDI ( var A : MATRIX ; lda, n : integer ; var IPvt : IVECTOR ; var Det : DVECTOR ; var Work : VECTOR ; JobCode : integer ) ; { } { PROGRAM: } { } { SGEDI - computes the determinant and inverse of } { a matrix using the factors computed by } { SGECO or SGEFA. } { } { VERSION: DATE: } { } { 1.0/TURBO Pascal 2.0 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 } { } { PASCAL translation; } { } { Adam Fritz } { 133 Main Street } { Afton, New York 13730 } { } const TEN : real = 10.0 ; iTEN : integer = 10 ; var d : real ; e : integer ; t : real ; i, j, k, l : integer ; kb, kp1, nm1 : integer ; begin { Compute Determinant } if (JobCode div iTEN) <> 0 then begin d := 1.0 ; e := 0 ; i := 0 ; while (i < n) and (d <> 0.0) do begin i := i + 1 ; if IPvt[i] <> i then d := -d ; d := A[i,i] * d ; if d <> 0.0 then if Abs(d) < 1.0 then repeat d := TEN * d ; e := e - 1 until Abs(d) >= 1.0 else if Abs(d) >= 10.0 then repeat d := d/TEN ; e := e + 1 until Abs(d) < 10.0 end ; Det[1] := d ; Det[2] := e end ; { Compute inverse(U) } if (JobCode mod iTEN) <> 0 then begin for k := 1 to n do begin A[k,k] := 1.0/A[k,k] ; t := -A[k,k] ; SSCAL(k-1, t, A[1,k], lda) ; if n > k then begin kp1 := k + 1 ; for j := kp1 to n do begin t := A[k,j] ; A[k,j] := 0.0 ; SAXPY(k, t, A[1,k], lda, A[1,j], lda) end end end ; { Form inverse(U) * inverse(L) } if n > 1 then begin nm1 := n - 1 ; for kb := 1 to nm1 do begin k := n - kb ; kp1 := k + 1 ; for i := kp1 to n do begin Work[i] := A[i,k] ; A[i,k] := 0.0 end ; for j := kp1 to n do begin t := Work[j] ; SAXPY(n, t, A[1,j], lda, A[1,k], lda) end ; l := IPvt[k] ; if l <> k then SSWAP(n, A[1,k], lda, A[1,l], lda) end end end end ; {Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730} {Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730} program main ; { } { Program: LINPACK SGEFA and SGESL Test Driver. } { } { Version: 1.0/T Date: 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 } { } {$I LINPACK.con CONSTANT Declarations } {$I LINPACK.typ TYPE Declarations } {$I LINPACK.var VARIABLE Declarations } {$I Hilgen.pas Test System Generator } {$I BLAS.pas Basic Linear Algebra } {$I SGEFA.pas LINPACK Factor } {$I SGESL.pas LINPACK Solve } {$I OUT.pas SICE Output Routine } begin { Initialize } writeln('LINPACK SGEFA and SGESL Test Program.') ; writeln ; { Get Order } n := 0 ; while (n < 1) or (n > lda) do begin write('Order: ') ; readln(n) end ; { Get Print Code } write('Print Code: ') ; readln (PrintCode) ; { Generate Test System } SYSGEN(A, lda, n, b) ; if PrintCode > 0 then begin writeln ; writeln('Original System (by column):') ; writeln ; OUT(A[1,1], lda, n, n) ; OUT(b[1], lda, n, 1) end ; { Factor the Matrix } writeln('begin: SGEFA') ; SGEFA(A, lda, n, IPvt, InfoCode) ; writeln('end: SGEFA') ; if PrintCode > 0 then begin writeln ; writeln('Factored Matrix (by column):') ; writeln ; OUT(A[1,1], lda, n, n) end ; { Solve the Linear System } if InfoCode = 0 then begin JobCode := 0 ; writeln('begin: SGESL') ; SGESL(A, lda, n, IPvt, b, JobCode) ; writeln('end: SGESL') ; if PrintCode > 0 then begin writeln ; writeln('Solution Vector:') ; writeln ; OUT(b[1], lda, n, 1) end end { Zero Diagonal Element Produced } else writeln('Error: Possible Singular System.') ; { Done } writeln('End of Test.') end. {Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730} {Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730} program main ; { } { Program: LINPACK SGECO and SGEFA Test Driver. } { } { Version: 1.0/T Date: 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 } { } {$I LINPACK.con CONSTANT Declarations } {$I LINPACK.typ TYPE Declarations } {$I LINPACK.var VARIABLE Declarations } {$I Hilgen.pas Test System Generator } {$I BLAS.pas Basic Linear Algebra } {$I SGEFA.pas LINPACK Factor } {$I SGECO.pas LINPACK Condition } {$I OUT.pas SICE Output Routine } begin { Initialize } writeln('LINPACK SGECO and SGEFA Test Program.') ; writeln ; { Get Order } n := 0 ; while (n < 1) or (n > lda) do begin write('Order: ') ; readln(n) end ; { Get Print Code } write('Print Code: ') ; readln (PrintCode) ; { Generate Test System } SYSGEN(A, lda, n, b) ; if PrintCode > 0 then begin writeln ; writeln('Original System (by column):') ; writeln ; OUT(A[1,1], lda, n, n) ; OUT(b[1], lda, n, 1) end ; { Compute the Condition } SGECO(A, lda, n, IPvt, RCond, Work) ; if PrintCode > 0 then begin writeln ; write('RCond: ',RCond:14) ; RCond := (1.0 + RCond) - 1.0 ; writeln(', RCond: ',RCond:14) ; writeln end ; { Done } writeln('End of Test.') end. {Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730} {Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730} program main ; { } { Program: LINPACK SGECO and SGEDI Test Driver. } { } { Version: 1.0/T Date: 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 } { } {$I LINPACK.con CONSTANT Declarations } {$I LINPACK.typ TYPE Declarations } {$I LINPACK.var VARIABLE Declarations } i, j : integer ; AA : MATRIX ; P : MATRIX ; {$I Hilgen.pas Test System Generator } {$I BLAS.pas Basic Linear Algebra } {$I SGEFA.pas LINPACK Factor } {$I SGECO.pas LINPACK Condition } {$I SGEDI.pas LINPACK Determinant and } { Inverse } {$I OUT.pas SICE Output Routine } begin { Initialize } writeln('LINPACK SGECO and SDEDI Test Program.') ; writeln ; { Get Order } n := 0 ; while (n < 1) or (n > lda) do begin write('Order: ') ; readln(n) end ; { Get Print Code } write('Print Code: ') ; readln (PrintCode) ; { Generate Test System } SYSGEN(A, lda, n, b) ; if PrintCode > 0 then begin writeln ; writeln('Original System (by column):') ; writeln ; OUT(A[1,1], lda, n, n) ; OUT(b[1], lda, n, 1) end ; { Copy Original by Rows } for i := 1 to n do SCOPY(n, A[i,1], 1, AA[i,1], 1) ; { Compute Condition } SGECO(A, lda, n, IPvt, RCond, Work) ; if PrintCode > 0 then begin writeln ; write('RCond: ',RCond:14) ; writeln(', RCond: ',((1.0+RCond)-1.0):14) ; writeln end ; { Compute Determinant and Inverse } if RCond <> 0.0 then begin JobCode := 11 ; SGEDI(A, lda, n, IPvt, Det, Work, JobCode) ; if PrintCode > 0 then begin writeln ; writeln('Inverse (by column):') ; writeln ; OUT(A[1,1], lda, n, n) ; writeln('Determinant: ',Det[1]:14:8,'E',Det[2]:3:0) end ; { Original Times Inverse } for i := 1 to n do for j := 1 to n do P[i,j] := SDOT(n, AA[i,1], 1, A[1,j], lda) ; if PrintCode > 0 then begin writeln ; writeln('Original times Inverse (by column):') ; writeln ; OUT(P[1,1], lda, n, n) end ; { Try to Restore Original } SGEFA(A, lda, n, IPvt, InfoCode) ; if InfoCode = 0 then begin SGEDI(A, lda, n, IPvt, Det, Work, JobCode) ; if PrintCode > 0 then begin writeln ; writeln('Inverse of Inverse (by column):') ; writeln ; OUT(A[1,1], lda, n, n) end end else writeln('Error: Inverse is Singular.') end { Zero Diagonal Element Produced } else writeln('Error: Original is Singular.') ; { Done } writeln('End of Test.') end. {Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730}