package linear_algebra;

import linear_algebra.*;
import corejava.Console;

import java.util.*;
import java.lang.*;
import java.io.*;

/**
*
*This class tests the (LINPACK) SVDC method.
*
*It
*
    *
  1. Randomly fills an n by n matrix, X. *
  2. Randomly generates a vector b0. *
  3. Calculates the vector y = X(b0). *
  4. Performs a UDV´ decomposition of X. *
  5. Solves the system Xb = y in an effort to recover b0. *
  6. Handles the regression equation y = a0 + a1*x + a2*x^2 + e * where a0 = a1 = a2 = 1. Its estimates of the parameters * should be close to 1. Its estimate of the standard deviation * should be close to the input value. *
  7. It tests the ability of the method to obtain the eigenvalues * and eigenvectors of a symmetric, positive definite matrix. *
* *@author Steve Verrill *@version .6 --- October 12, 1997 * */ public class SVDCTest_f77 extends Object { public static void main (String args[]) { /* Some of the variables of main: */ long randstart; int n,p; int i,j,k; double randlow,randhigh,randdiff; double x[][] = new double[101][101]; double copyx[][] = new double[101][101]; double y[] = new double[101]; double b0[] = new double[101]; double b[] = new double[101]; double s[] = new double[101]; double e[] = new double[101]; double u[][] = new double[101][101]; double ut[][] = new double[101][101]; double v[][] = new double[101][101]; double work[] = new double[101]; double pred[] = new double[101]; double uty[] = new double[101]; double sd,rss,diff; double onedsqrt4; double diag[][] = new double[5][5]; double dut[][] = new double[5][5]; int job; /* Console is a public domain class described in Cornell and Horstmann's Core Java (SunSoft Press, Prentice-Hall). */ n = Console.readInt("\nWhat is the n value? (100 or less) "); randstart = Console.readInt("\nWhat is the starting value (an " + "integer)\n" + "for the random number generator? "); Random rand01 = new Random(randstart); randlow = Console.readDouble("\nWhat is randlow? "); randhigh = Console.readDouble("\nWhat is randhigh? "); randdiff = randhigh - randlow; sd = Console.readDouble("\nWhat is the standard deviation? "); /* Generate an n by n X */ for (i = 1; i <= n; i++) { for (j = 1; j <= n; j++) { x[i][j] = randlow + rand01.nextDouble()*randdiff; } } /* System.out.print("\n\n"); System.out.print("*********************************************\n"); System.out.print("\nThe X matrix is \n\n"); for (i = 1; i <= n; i++) { for (j = 1; j <= n; j++) { System.out.print(x[i][j] + " "); } System.out.print("\n"); } */ /* Generate a vector b */ for (i = 1; i <= n; i++) { b0[i] = randlow + rand01.nextDouble()*randdiff; } /* Calculate Xb = y */ for (i = 1; i <= n; i++) { y[i] = 0.0; for (j = 1; j <= n; j++) { y[i] += x[i][j]*b0[j]; } } job = 11; try { SVDC_f77.dsvdc_f77(x,n,n,s,e,u,v,work,job); /* This should obtain the solution to Xb = y. */ Blas_f77.mattran_f77(u,ut,n,n); Blas_f77.matvec_f77(ut,y,uty,n,n); for (i = 1; i <= n; i++) { uty[i] = uty[i]/s[i]; } Blas_f77.matvec_f77(v,uty,b,n,n); System.out.print("\nThe input b vector was:\n\n"); for (i = 1; i <= n; i++) { System.out.print(b0[i] + "\n"); } System.out.print("\n\nThe recovered b vector was:\n\n"); for (i = 1; i <= n; i++) { System.out.print(b[i] + "\n"); } } catch (SVDCException svdce) { System.out.print("\nThere was an SVDC_f77 exception on the" + " first test.\n\n" + "The info value from SVDC_f77 was " + svdce.info + ".\n\n"); } System.out.print("\n"); System.out.print("\n"); System.out.print("*********************************************\n"); System.out.print("\n"); /* Generate a new X matrix and the y vector */ for (i = 1; i <= n; i++) { x[i][1] = 1.0; copyx[i][1] = x[i][1]; x[i][2] = -1.0 + rand01.nextDouble()*2.0; copyx[i][2] = x[i][2]; x[i][3] = x[i][2]*x[i][2]; copyx[i][3] = x[i][3]; y[i] = x[i][1] + x[i][2] + x[i][3] + sd*normi(rand01.nextDouble()); } job = 21; try { SVDC_f77.dsvdc_f77(x,n,3,s,e,u,v,work,job); /* This should obtain the solution to Xb = y. */ Blas_f77.mattran_f77(u,ut,n,3); Blas_f77.matvec_f77(ut,y,uty,3,n); for (i = 1; i <= 3; i++) { uty[i] = uty[i]/s[i]; } b[1] = 0.0; b[2] = 0.0; b[3] = 0.0; Blas_f77.matvec_f77(v,uty,b,3,3); System.out.print("\nThe input b vector was: \n\n1\n1\n1\n\n"); System.out.print("\n\nThe recovered b vector was:\n\n"); for (i = 1; i <= 3; i++) { System.out.print(b[i] + "\n"); } Blas_f77.matvec_f77(copyx,b,pred,n,3); rss = 0.0; for (i = 1; i <= n; i++) { diff = y[i] - pred[i]; rss += diff*diff; } sd = Math.sqrt(rss/(n - 1.0)); System.out.print("\nThe estimate of sigma from the residuals"); System.out.print(" was " + sd + "\n"); } catch (SVDCException svdce) { System.out.print("\nThere was an SVDC_f77 exception on the" + " second test.\n\n" + "The info value from SVDC_f77 was " + svdce.info + ".\n\n"); } System.out.print("\n"); System.out.print("\n"); System.out.print("*********************************************\n"); System.out.print("\n\nEigenvalue, eigenvector test\n\n"); /* Generate a symmetric, positive definite matrix with eigenvalues 4,3,2,1. */ for (i = 1; i <= 4; i++) { y[i] = 1.0; } onedsqrt4 = 1.0/Blas_f77.dnrm2_f77(4,y,1); for (i = 1; i <= 4; i++) { u[i][1] = onedsqrt4; } u[1][2] = onedsqrt4; u[2][2] = -onedsqrt4; u[3][2] = onedsqrt4; u[4][2] = -onedsqrt4; u[1][3] = onedsqrt4; u[2][3] = -onedsqrt4; u[3][3] = -onedsqrt4; u[4][3] = onedsqrt4; u[1][4] = onedsqrt4; u[2][4] = onedsqrt4; u[3][4] = -onedsqrt4; u[4][4] = -onedsqrt4; Blas_f77.mattran_f77(u,ut,4,4); for (i = 1; i <= 4; i++) { for (j = 1; j <= 4; j++) { diag[i][j] = 0.0; } } diag[1][1] = 4.0; diag[2][2] = 3.0; diag[3][3] = 2.0; diag[4][4] = 1.0; Blas_f77.matmat_f77(diag,ut,dut,4,4,4); Blas_f77.matmat_f77(u,dut,x,4,4,4); job = 10; try { SVDC_f77.dsvdc_f77(x,4,4,s,e,u,v,work,job); System.out.print("\nThe eigenvalues should be 4,3,2,1.\n"); System.out.print("\nThe calculated eigenvalues are:\n\n"); System.out.print(s[1] + "\n"); System.out.print(s[2] + "\n"); System.out.print(s[3] + "\n"); System.out.print(s[4] + "\n"); System.out.print("\n\nThe eigenvectors should be (up to sign):\n\n"); System.out.print(".5 .5 .5 .5\n"); System.out.print(".5 -.5 -.5 .5\n"); System.out.print(".5 .5 -.5 -.5\n"); System.out.print(".5 -.5 .5 -.5\n"); System.out.print("\n\nThe calculated eigenvectors are:\n\n"); System.out.print(u[1][1] + " " + u[1][2] + " " + u[1][3] + " " + u[1][4] + "\n"); System.out.print(u[2][1] + " " + u[2][2] + " " + u[2][3] + " " + u[2][4] + "\n"); System.out.print(u[3][1] + " " + u[3][2] + " " + u[3][3] + " " + u[3][4] + "\n"); System.out.print(u[4][1] + " " + u[4][2] + " " + u[4][3] + " " + u[4][4] + "\n"); } catch (SVDCException svdce) { System.out.print("\nThere was an SVDC_f77 exception on the" + " eigenvalue, eigenvector test.\n\n" + "The info value from SVDC_f77 was " + svdce.info + ".\n\n"); } } /** * *

*This is a normal cdf inverse routine. * *Created by Steve Verrill, March 1997. * *@param u The value (between 0 and 1) to invert. * */ public static double normi (double u) { // This is a temporary and potentially ugly way to generate // normal random numbers from uniform random numbers double c[] = new double[4]; double d[] = new double[4]; double normi,arg,t,t2,t3,xnum,xden; c[1] = 2.515517; c[2] = .802853; c[3] = .010328; d[1] = 1.432788; d[2] = .189269; d[3] = .001308; if (u <= .5) { arg = -2.0*Math.log(u); t = Math.sqrt(arg); t2 = t*t; t3 = t2*t; xnum = c[1] + c[2]*t + c[3]*t2; xden = 1.0 + d[1]*t + d[2]*t2 + d[3]*t3; normi = -(t - xnum/xden); return normi; } else { arg = -2.0*Math.log(1.0 - u); t = Math.sqrt(arg); t2 = t*t; t3 = t2*t; xnum = c[1] + c[2]*t + c[3]*t2; xden = 1.0 + d[1]*t + d[2]*t2 + d[3]*t3; normi = t - xnum/xden; return normi; } } }