package linear_algebra;

import java.util.*;
import java.lang.*;
import java.io.*;
import linear_algebra.*;
import corejava.Console;

/**
*
*This class tests the LU_f77 methods.
*
*It
*
    *
  1. Randomly generates numbers between randlow and randhigh *and fills a matrix A with them. *
  2. Randomly generates a vector x. *
  3. Calculates the vector b = Ax. *
  4. Performs a LU decomposition of A. *(This tests LU_f77.dgefa_f77.) *
  5. Solves the system Az = b in an effort to recover x. *(This tests LU_f77.dgesl_f77.) *
  6. Obtains an estimate of A^{-1} and compares AA^{-1} with the *identity matrix. (This tests the matrix inverse portion *of LU_f77.dgedi_f77.) *
  7. Via LU_f77.dgefa_f77 and LU_f77.dgedi_f77, obtains the *determinant, detval2, of transpose(R)R where R is upper triangular, *and compares it with detval1 = the square of the product of the *diagonal elements of R. (This tests the determinant portion *of LU_f77.dgedi_f77.) *
* *@author Steve Verrill *@version .5 --- March 11, 1998 * */ public class LUTest_f77 extends Object { public static void main (String args[]) { /* Some of the variables of main: 1.) randstart --- The integer starting value for the random number generator. 2.) n --- The order of the vectors and matrices. 3.) a[][] --- Initially, it holds all of A. After dgefa is applied, it contains "an upper triangular matrix and the multipliers which were used to obtain it." 4.) a2[][] --- A copy of A. 5.) r[][] --- Used in test of determinant calculation. 6.) rt[][] --- Used in test of determinant calculation. 7.) detval1 --- Used in test of determinant calculation. 8.) detval2 --- Used in test of determinant calculation. 9.) x[] --- The x in Ax = b. 10.) work[] --- Work array needed to help obtain A^[-1]. 11.) b[] --- The b in Ax = b. 12.) ipvt[] --- a vector of pivot indices 13.) iden[][] --- should be the identity matrix after AA^[-1] is calculated */ long randstart; int n; int i,j,k; double randlow,randhigh,randdiff; double a[][] = new double[101][101]; double a2[][] = new double[101][101]; double x[] = new double[101]; double work[] = new double[101]; double b[] = new double[101]; double det[] = new double[3]; int ipvt[] = new int[101]; double iden[][] = new double[101][101]; double r[][] = new double[101][101]; double rt[][] = new double[101][101]; double detval1,detval2; /* Console is a public domain class described in Cornell and Horstmann's Core Java (SunSoft Press, Prentice-Hall). */ randstart = Console.readInt("\nWhat is the starting value (an " + "integer)\n" + "for the random number generator? "); n = Console.readInt("\nWhat is the n value? (100 or less) "); Random rand01 = new Random(randstart); randlow = Console.readDouble("\nWhat is randlow? "); randhigh = Console.readDouble("\nWhat is randhigh? "); randdiff = randhigh - randlow; /* Generate a matrix. */ for (i = 1; i <= n; i++) { for (j = 1; j <= n; j++) { a[i][j] = randlow + rand01.nextDouble()*randdiff; a2[i][j] = a[i][j]; r[i][j] = 0.0; } } /* Material needed to test determinant calculations */ detval1 = 1.0; for (i = 1; i <= n; i++) { for (j = i; j <= n; j++) { r[i][j] = a[i][j]; } detval1 *= r[i][i]; } detval1 = detval1*detval1; Blas_f77.mattran_f77(r,rt,n,n); /* System.out.print("\n\n"); System.out.print("*********************************************\n"); System.out.print("\nThe a matrix is \n\n"); for (i = 1; i <= n; i++) { for (j = 1; j <= n; j++) { System.out.print(a[i][j] + " "); } System.out.print("\n"); } */ /* Generate a vector x */ for (i = 1; i <= n; i++) { x[i] = randlow + rand01.nextDouble()*randdiff; } /* Calculate Ax = b */ for (i = 1; i <= n; i++) { b[i] = 0.0; for (j = 1; j <= n; j++) { b[i] += a[i][j]*x[j]; } } /* Factor A. */ try { LU_f77.dgefa_f77(a,n,ipvt); } catch (NotFullRankException nfr) { System.out.print("\nThe matrix A was not of full rank" + " so it was not possible\n" + "to test the LU decomposition.\n\n"); System.out.print("\nThe info value was " + nfr.info + ".\n\n"); System.exit(0); } System.out.print("\n"); System.out.print("*********************************************\n"); System.out.print("\n"); /* This should obtain the solution to Ax = b. The solution x will be in the vector b. */ LU_f77.dgesl_f77(a,n,ipvt,b,0); System.out.print("\nThe input x vector was\n\n"); for (i = 1; i <= n; i++) { System.out.print(x[i] + "\n"); } System.out.print("\n\nThe recovered x vector was\n\n"); for (i = 1; i <= n; i++) { System.out.print(b[i] + "\n"); } System.out.print("\n"); System.out.print("\n"); System.out.print("*********************************************\n"); System.out.print("\n"); /* Obtain the inverse of A. */ LU_f77.dgedi_f77(a,n,ipvt,det,work,1); /* now a*a2 should equal I */ for (i = 1; i <= n; i++) { for (j = 1; j <= n; j++) { iden[i][j] = 0.0; for (k = 1; k <= n; k++) { iden[i][j] += a2[i][k]*a[k][j]; } } } System.out.print("\nA test of LU_f77.dgedi_f77's calculation" + " of inverses\n\n" + "The I matrix is:\n\n"); for (i = 1; i <= n; i++) { for (j = 1; j <= n; j++) { System.out.print(iden[i][j] + " "); } System.out.print("\n"); } System.out.print("\n"); System.out.print("\n"); System.out.print("*********************************************\n"); System.out.print("\n"); /* Obtain the determinant of (R^T)(R). */ Blas_f77.matmat_f77(rt,r,a,n,n,n); try { LU_f77.dgefa_f77(a,n,ipvt); } catch (NotFullRankException nfr) { System.out.print("\nThe matrix RT*R was not of full rank" + " so the determinant should be 0.\n"); System.out.print("\nThe info value was " + nfr.info + ".\n\n"); } LU_f77.dgedi_f77(a,n,ipvt,det,work,10); detval2 = det[1]*Math.pow(10.0,det[2]); System.out.print("\nA test of LU_f77.dgedi_f77's calculation" + " of determinants\n\n"); System.out.print("detval1 = " + detval1 + "\n"); System.out.print("detval2 = " + detval2 + "\n"); System.out.print("\n"); System.out.print("\n"); System.out.print("*********************************************\n"); System.out.print("\n"); } }