```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
*
* Randomly fills an n by n matrix, X.
* Randomly generates a vector b0.
* Calculates the vector y = X(b0).
* Performs a UDV´ decomposition of X.
* Solves the system Xb = y in an effort to recover b0.
* 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.
* 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;
double copyx[][] = new double;

double y[] = new double;
double b0[] = new double;
double b[] = new double;

double s[] = new double;
double e[] = new double;
double u[][] = new double;
double ut[][] = new double;
double v[][] = new double;
double work[] = new double;

double pred[] = new double;

double uty[] = new double;

double onedsqrt4;

double diag[][] = new double;
double dut[][] = new double;

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.0;
copyx[i] = x[i];
x[i] = -1.0 + rand01.nextDouble()*2.0;
copyx[i] = x[i];
x[i] = x[i]*x[i];
copyx[i] = x[i];

y[i] = x[i] + x[i] + x[i] + 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 = 0.0;
b = 0.0;
b = 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);

for (i = 1; i <= n; i++) {

diff = y[i] - pred[i];

}

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] = onedsqrt4;

}

u = onedsqrt4;
u = -onedsqrt4;
u = onedsqrt4;
u = -onedsqrt4;

u = onedsqrt4;
u = -onedsqrt4;
u = -onedsqrt4;
u = onedsqrt4;

u = onedsqrt4;
u = onedsqrt4;
u = -onedsqrt4;
u = -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 = 4.0;
diag = 3.0;
diag = 2.0;
diag = 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 + "\n");
System.out.print(s + "\n");
System.out.print(s + "\n");
System.out.print(s + "\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 + "  " + u + "  " +
u + "  " + u + "\n");
System.out.print(u + "  " + u + "  " +
u + "  " + u + "\n");
System.out.print(u + "  " + u + "  " +
u + "  " + u + "\n");
System.out.print(u + "  " + u + "  " +
u + "  " + u + "\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;
double d[] = new double;

double normi,arg,t,t2,t3,xnum,xden;

c = 2.515517;
c = .802853;
c = .010328;

d = 1.432788;
d = .189269;
d = .001308;

if (u <= .5) {

arg = -2.0*Math.log(u);
t = Math.sqrt(arg);
t2 = t*t;
t3 = t2*t;

xnum = c + c*t + c*t2;
xden = 1.0 + d*t + d*t2 + d*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 + c*t + c*t2;
xden = 1.0 + d*t + d*t2 + d*t3;

normi = t - xnum/xden;

return normi;

}

}

}
```