/*
    LU_f77.java copyright claim:

    This software is based on public domain LINPACK routines.
    It was translated from FORTRAN to Java by a US government employee 
    on official time.  Thus this software is also in the public domain.


    The translator's mail address is:

    Steve Verrill 
    USDA Forest Products Laboratory
    1 Gifford Pinchot Drive
    Madison, Wisconsin
    53705


    The translator's e-mail address is:

    steve@swst.org


***********************************************************************

DISCLAIMER OF WARRANTIES:

THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND. 
THE TRANSLATOR DOES NOT WARRANT, GUARANTEE OR MAKE ANY REPRESENTATIONS 
REGARDING THE SOFTWARE OR DOCUMENTATION IN TERMS OF THEIR CORRECTNESS, 
RELIABILITY, CURRENTNESS, OR OTHERWISE. THE ENTIRE RISK AS TO 
THE RESULTS AND PERFORMANCE OF THE SOFTWARE IS ASSUMED BY YOU. 
IN NO CASE WILL ANY PARTY INVOLVED WITH THE CREATION OR DISTRIBUTION 
OF THE SOFTWARE BE LIABLE FOR ANY DAMAGE THAT MAY RESULT FROM THE USE 
OF THIS SOFTWARE.

Sorry about that.

***********************************************************************


History:

Date        Translator        Changes

3/11/98     Steve Verrill     Translated

*/


package linear_algebra;

import linear_algebra.*;


/**
*
*

*This class contains the LINPACK DGEFA (LU factorization), *DGESL (solve), and DGEDI (determinant and inverse) routines. * *

*IMPORTANT: The "_f77" suffixes indicate that these routines use *FORTRAN style indexing. For example, you will see * * for (i = 1; i <= n; i++) * *rather than * * for (i = 0; i < n; i++) * *To use the "_f77" routines you will have to declare your vectors *and matrices to be one element larger (e.g., v[101] rather than *v[100], and a[101][101] rather than a[100][100]), and you will have *to fill elements 1 through n rather than elements 0 through n - 1. *Versions of these programs that use C/Java style indexing are *also available. They end with the suffix "_j". * *

*This class was translated by a statistician from FORTRAN versions of *the LINPACK routines. It is NOT an official translation. When *public domain Java numerical analysis routines become available *from the people who produce LAPACK, then THE CODE PRODUCED *BY THE NUMERICAL ANALYSTS SHOULD BE USED. * *

*Meanwhile, if you have suggestions for improving this *code, please contact Steve Verrill at sverrill@fs.fed.us. * *@author Steve Verrill *@version .5 --- March, 1998 * */ public class LU_f77 extends Object { /** * *

*This method decomposes an n by n *matrix A into a product, LU, where L is a lower triangular matrix *and U is an upper triangular matrix. *For details, see the comments in the code. *This method is a translation from FORTRAN to Java of the LINPACK subroutine *DGEFA. In the LINPACK listing DGEFA is attributed to Cleve Moler *with a date of 8/14/78. * *Translated by Steve Verrill, March 10, 1998. * *@param a The matrix to be decomposed *@param n The order of a *@param ipvt A vector of pivot indices * */ public static void dgefa_f77 (double a[][], int n, int ipvt[]) throws NotFullRankException { /* Here is a copy of the LINPACK documentation (from the SLATEC version): C***BEGIN PROLOGUE DGEFA C***DATE WRITTEN 780814 (YYMMDD) C***REVISION DATE 861211 (YYMMDD) C***CATEGORY NO. D2A1 C***KEYWORDS LIBRARY=SLATEC(LINPACK), C TYPE=DOUBLE PRECISION(SGEFA-S DGEFA-D CGEFA-C), C GENERAL MATRIX,LINEAR ALGEBRA,MATRIX,MATRIX FACTORIZATION C***AUTHOR MOLER, C. B., (U. OF NEW MEXICO) C***PURPOSE Factors a double precision matrix by Gaussian elimination. C***DESCRIPTION C C DGEFA factors a double precision matrix by Gaussian elimination. C C DGEFA is usually called by DGECO, but it can be called C directly with a saving in time if RCOND is not needed. C (Time for DGECO) = (1 + 9/N)*(Time for DGEFA) . C C On Entry C C A DOUBLE PRECISION(LDA, N) C the matrix to be factored. C C LDA INTEGER C the leading dimension of the array A . C C N INTEGER C the order of the matrix A . C C On Return C C A an upper triangular matrix and the multipliers C which were used to obtain it. C The factorization can be written A = L*U where C L is a product of permutation and unit lower C triangular matrices and U is upper triangular. C C IPVT INTEGER(N) C an integer vector of pivot indices. C C INFO INTEGER C = 0 normal value. C = K if U(K,K) .EQ. 0.0 . This is not an error C condition for this subroutine, but it does C indicate that DGESL or DGEDI will divide by zero C if called. Use RCOND in DGECO for a reliable C indication of singularity. For the Java version, info is returned as an argument to NotFullRankException. C C LINPACK. This version dated 08/14/78 . C Cleve Moler, University of New Mexico, Argonne National Lab. C C Subroutines and Functions C C BLAS DAXPY,DSCAL,IDAMAX C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W., C *LINPACK USERS GUIDE*, SIAM, 1979. C***ROUTINES CALLED DAXPY,DSCAL,IDAMAX C***END PROLOGUE DGEFA */ double t; int isamax,j,k,kp1,l,nm1; // Gaussian elimination with partial pivoting nm1 = n - 1; // if (nm1 >= 1) { for (k = 1; k <= nm1; k++) { kp1 = k + 1; // Find l = pivot index l = Blas_f77.colisamax_f77(n-k+1,a,1,k,k) + k - 1; ipvt[k] = l; // Zero pivot implies this column already 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]; Blas_f77.colscal_f77(n-k,t,a,k+1,k); // 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; } Blas_f77.colaxpy_f77(n-k,t,a,k+1,k,j); } } else { throw new NotFullRankException(k); } } // } ipvt[n] = n; if (a[n][n] == 0.0) { throw new NotFullRankException(n); } return; } /** * *

*This method uses the LU decomposition provided by *DGEFA to solve the equation Ax = b where A is a full *rank n by n matrix. *For details, see the comments in the code. *This method is a translation from FORTRAN to Java of the LINPACK subroutine *DGESL. In the LINPACK listing DGESL is attributed to Cleve Moler *with a date of 8/14/78. * *Translated by Steve Verrill, March 11, 1998. * *@param a a[][] *@param n The order of a *@param ipvt A vector of pivot indices *@param b Input --- the vector b in Ax = b, Output --- the * vector x in Ax = b *@param job 0 --- solve Ax = b, nonzero --- solve Transpose(A)x = b * */ public static void dgesl_f77 (double a[][], int n, int ipvt[], double b[], int job) { /* Here is a copy of the LINPACK documentation (from the SLATEC version): C***BEGIN PROLOGUE DGESL C***DATE WRITTEN 780814 (YYMMDD) C***REVISION DATE 861211 (YYMMDD) C***CATEGORY NO. D2A1 C***KEYWORDS LIBRARY=SLATEC(LINPACK), C TYPE=DOUBLE PRECISION(SGESL-S DGESL-D CGESL-C), C LINEAR ALGEBRA,MATRIX,SOLVE C***AUTHOR MOLER, C. B., (U. OF NEW MEXICO) C***PURPOSE Solves the double precision system A*X=B or TRANS(A)*X=B C using the factors computed by DGECO or DGEFA. C***DESCRIPTION C C DGESL solves the double precision system C A * X = B or TRANS(A) * X = B C using the factors computed by DGECO or DGEFA. C C On Entry C C A DOUBLE PRECISION(LDA, N) C the output from DGECO or DGEFA. C C LDA INTEGER C the leading dimension of the array A . C C N INTEGER C the order of the matrix A . C C IPVT INTEGER(N) C the pivot vector from DGECO or DGEFA. C C B DOUBLE PRECISION(N) C the right hand side vector. C C JOB INTEGER C = 0 to solve A*X = B , C = nonzero to solve TRANS(A)*X = B where C TRANS(A) is the transpose. C C On Return C C B the solution vector X . C C Error Condition C C A division by zero will occur if the input factor contains a C zero on the diagonal. Technically this indicates singularity C but it is often caused by improper arguments or improper C setting of LDA . It will not occur if the subroutines are C called correctly and if DGECO has set RCOND .GT. 0.0 C or DGEFA has set INFO .EQ. 0 . C C To compute INVERSE(A) * C where C is a matrix C with P columns C CALL DGECO(A,LDA,N,IPVT,RCOND,Z) C IF (RCOND is too small) GO TO ... C DO 10 J = 1, P C CALL DGESL(A,LDA,N,IPVT,C(1,J),0) C 10 CONTINUE C C LINPACK. This version dated 08/14/78 . C Cleve Moler, University of New Mexico, Argonne National Lab. C C Subroutines and Functions C C BLAS DAXPY,DDOT C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W., C *LINPACK USERS GUIDE*, SIAM, 1979. C***ROUTINES CALLED DAXPY,DDOT C***END PROLOGUE DGESL */ double t; int k,kb,l,nm1; nm1 = n - 1; if (job == 0) { // job = 0, solve A*x = b // First solve L*y = b // if (nm1 >= 1) { for (k = 1; k <= nm1; k++) { l = ipvt[k]; t = b[l]; if (l != k) { b[l] = b[k]; b[k] = t; } Blas_f77.colvaxpy_f77(n-k,t,a,b,k+1,k); } // } // Now solve U*x = y for (kb = 1; kb <= n; kb++) { k = n + 1 - kb; b[k] /= a[k][k]; t = -b[k]; Blas_f77.colvaxpy_f77(k-1,t,a,b,1,k); } } else { // job = nonzero, solve transpose(A)*x = b // First solve transpose(U)*y = b for (k = 1; k <= n; k++) { t = Blas_f77.colvdot_f77(k-1,a,b,1,k); b[k] = (b[k] - t)/a[k][k]; } // Now solve transpose(L)*x = y // if (nm1 >= 1) { for (kb = 1; kb <= nm1; kb++) { k = n - kb; b[k] += Blas_f77.colvdot_f77(n-k,a,b,k+1,k); l = ipvt[k]; if (l != k) { t = b[l]; b[l] = b[k]; b[k] = t; } } // } } return; } /** * *

*This method uses the LU decomposition provided by *DGEFA to obtain the determinant and/or inverse of a *full rank n by n matrix. *For details, see the comments in the code. *This method is a translation from FORTRAN to Java of the LINPACK subroutine *DGEDI. In the LINPACK listing DGEDI is attributed to Cleve Moler *with a date of 8/14/78. * *Translated by Steve Verrill, March 11, 1998. * *@param a a[][] *@param n The order of a *@param ipvt A vector of pivot indices *@param det det[] *@param work work[] *@param job Indicates whether a determinant, inverse, * or both is desired * */ public static void dgedi_f77 (double a[][], int n, int ipvt[], double det[], double work[], int job) { /* Here is a copy of the LINPACK documentation (from the SLATEC version): C***BEGIN PROLOGUE DGEDI C***DATE WRITTEN 780814 (YYMMDD) C***REVISION DATE 861211 (YYMMDD) C***CATEGORY NO. D3A1,D2A1 C***KEYWORDS LIBRARY=SLATEC(LINPACK), C TYPE=DOUBLE PRECISION(SGEDI-S DGEDI-D CGEDI-C), C DETERMINANT,INVERSE,LINEAR ALGEBRA,MATRIX, C MATRIX FACTORIZATION C***AUTHOR MOLER, C. B., (U. OF NEW MEXICO) C***PURPOSE Computes the determinant and inverse of a matrix using C factors computed by DGECO or DGEFA. C***DESCRIPTION C C DGEDI computes the determinant and inverse of a matrix C using the factors computed by DGECO or DGEFA. C C On Entry C C A DOUBLE PRECISION(LDA, N) C the output from DGECO or DGEFA. C C LDA INTEGER C the leading dimension of the array A . C C N INTEGER C the order of the matrix A . C C IPVT INTEGER(N) C the pivot vector from DGECO or DGEFA. C C WORK DOUBLE PRECISION(N) C work vector. Contents destroyed. C C JOB INTEGER C = 11 both determinant and inverse. C = 01 inverse only. C = 10 determinant only. C C On Return C C A inverse of original matrix if requested. C Otherwise unchanged. C C DET DOUBLE PRECISION(2) C determinant of original matrix if requested. C Otherwise not referenced. C Determinant = DET(1) * 10.0**DET(2) C with 1.0 .LE. DABS(DET(1)) .LT. 10.0 C or DET(1) .EQ. 0.0 . C C Error Condition C C A division by zero will occur if the input factor contains C a zero on the diagonal and the inverse is requested. C It will not occur if the subroutines are called correctly C and if DGECO has set RCOND .GT. 0.0 or DGEFA has set C INFO .EQ. 0 . C C LINPACK. This version dated 08/14/78 . C Cleve Moler, University of New Mexico, Argonne National Lab. C C Subroutines and Functions C C BLAS DAXPY,DSCAL,DSWAP C Fortran DABS,MOD C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W., C *LINPACK USERS GUIDE*, SIAM, 1979. C***ROUTINES CALLED DAXPY,DSCAL,DSWAP C***END PROLOGUE DGEDI */ double t,ten; int i,j,k,kb,kp1,l,nm1; if (job/10 != 0) { // Compute determinant det[1] = 1.0; det[2] = 0.0; ten = 10.0; for (i = 1; i <= n; i++) { if (ipvt[i] != i) det[1] = -det[1]; det[1] *= a[i][i]; if (det[1] == 0.0) break; while (Math.abs(det[1]) < 1.0) { det[1] *= ten; det[2]--; } while (Math.abs(det[1]) >= ten) { det[1] /= ten; det[2]++; } } } if ((job % 10) != 0) { // Compute inverse for (k = 1; k <= n; k++) { a[k][k] = 1.0/a[k][k]; t = -a[k][k]; Blas_f77.colscal_f77(k-1,t,a,1,k); kp1 = k + 1; if (n >= kp1) { for (j = kp1; j <= n; j++) { t = a[k][j]; a[k][j] = 0.0; Blas_f77.colaxpy_f77(k,t,a,1,k,j); } } } // Form inverse(U)*inverse(L) nm1 = n - 1; // if (nm1 >= 1) { for (kb = 1; kb <= nm1; kb++) { k = n - kb; 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]; Blas_f77.colaxpy_f77(n,t,a,1,j,k); } l = ipvt[k]; if (l != k) Blas_f77.colswap_f77(n,a,k,l); } // } } return; } }