Cholesky Decomposition

[Eric Gamess <egamess@kanaima.ciens.ucv.ve>]

Hello,

I did a program for Cholesky decomposition with Jama, but a have some
problems. The program is easy. It generates a random matrix A and
a random vector solution called xRandom. Them, it computes the
right-hand-side vector called b (b <== A*xRandom). xComputed is
a vector that receives the computed
solution (xComputed <== A.chol().solve(b)). Then, I show the norm of
xRandom-xComputed and have a result far from 0.0.

My program is included in the mail. Please, could you look at it and tell
me if you see and error.

Thank you very much.

Eric.



//////////////////////////////////////////////////////////////////////
import Jama.*;
 
 
class Problem
{
 
   static public void main(String args[])
   {
   int size;
 
   EasyIn easy=new EasyIn();
 
   System.out.print("Enter matrix size ... ");
   size=easy.readInt();
 
   Matrix A=new Matrix(size, size);
   Matrix xRandom=new Matrix(size, 1);
 
   create_problem(A, xRandom);
   Matrix b=A.times(xRandom);
 
   Matrix xComputed=null;
   try
      {
      xComputed=A.chol().solve(b);
      // Resolves the symmetric positive definite system of
      // linear equations
      }
   catch(RuntimeException re)
      {
      System.err.println("Non positive definite matrix");
      System.exit(1);
      }
 
   Matrix diff=xComputed.minus(xRandom);
   // Subtracts xRamdom from xComputed and, puts the result in diff
 
   double norm=diff.norm2();
   // Computes norm of vector diff
 
   System.out.println("Norm of difference: "+norm);
   System.exit(0);
   }
 
 
   static void create_problem(Matrix A, Matrix xRamdom)
   {
   int size=A.getRowDimension();
 
   double ATmp[][]=A.getArray();
   for(int i=0; i<size; i++)
      for(int j=i; j<size; j++)
         ATmp[i][j]=ATmp[j][i]=Math.random()*2.0-1.0;
   // The matrix is symmetric
 
   for(int i=0; i<size; i++)
      ATmp[i][i]+=size;
   // To have a better chance that the matrix be positive-definite,
   // I add a big value to elements of the main diagonal

   double xRamdomTmp[][]=xRamdom.getArray();
   for(int i=0; i<size; i++)
      xRamdomTmp[i][0]=Math.random()*2.0-1.0;
   } 
}