Example Codes of dsymatrix LAPACK Functions

//=============================================================================
/*! dsysv_check */
void dsysv_check_vector()
{
  cout << "############ check dsysv vector ############" << endl;
  
  srand(unsigned(time(NULL)));
  int N(3);
  
  //// make dsymatrix A  ////
  CPPL::dsymatrix A(N);
  for(int i=0; i<A.n; i++){
    for(int j=0; j<A.n; j++){
      A(i,j) =double( rand() /(RAND_MAX/10) );
    }
  }
  
  //// make dcovector b ////
  CPPL::dcovector b(N);
  for(int i=0; i<b.l; i++){
    b(i) =double( rand() /(RAND_MAX/10) );
  }
  
  //// make A_original and b_original ////
  CPPL::dsymatrix A_original(A);
  CPPL::dcovector b_original(b);
  cout << "A_original=\n" << A_original << endl;
  cout << "b_original=\n" << b_original << endl;
  
  //// solve Ax=b ////
  A.dsysv(b);
  
  //// print A, b and A_original*y ////
  cout << "modified A=\n" << A << endl;
  cout << "x=\n" << b << endl;
  cout << "b_original -A_original*x = (should be zero)\n" << b_original -A_original*b << endl;
}

//=============================================================================
void dsysv_check_matrix()
{
  cout << "############ check dsysv matrix ############" << endl;
  
  srand(unsigned(time(NULL)));
  int N(3);
  

  //// make dsymatrix A  ////
  CPPL::dsymatrix A(N);
  for(int i=0; i<A.n; i++){
    for(int j=0; j<A.n; j++){
      A(i,j) =double( rand() /(RAND_MAX/10) );
    }
  }
   
  //// make dsymatrix B  ////
  CPPL::dgematrix B(N,N);
  for(int i=0; i<B.m; i++){
    for(int j=0; j<B.n; j++){
      B(i,j) =double( rand() /(RAND_MAX/10) );
    }
  }
   
  //// make A_original and B_original ////
  CPPL::dsymatrix A_original(A);
  CPPL::dgematrix B_original(B);
  cout << "A_original=\n" << A_original << endl;
  cout << "B_original=\n" << B_original << endl;
  
  //// solve A*X=B ////
  A.dsysv(B);
  
  //// print A, B and A_original*B ////
  cout << "modified A=\n" << A << endl;
  cout << "X=\n" << B << endl;
  cout << "B_original -A_original*X = (should be zero)\n" << B_original -A_original*B << endl;
}

//=============================================================================
void dsyev_check_value()
{
  cout << "############ check dsyev value ############" << endl;
  
  srand(unsigned(time(NULL)));
  int N(3);
  
  //// make dsymatrix A  ////
  CPPL::dsymatrix A(N);
  for(int i=0; i<N; i++){ for(int j=0; j<N; j++){
    A(i,j) =double( rand() /(RAND_MAX/10) );
  }}
  
  //// make wr wi vr ////
  vector<double> w;
  
  //// make A_original ////
  CPPL::dsymatrix A_original(A);
  
  //// dsyev ////
  A.dsyev(w);
  //A.dsyev(w,1);
  
  //// print ////
  cout << "A=\n" << A << endl;
  cout << "A_original=\n" << A_original << endl;
  for(int i=0; i<N; i++){
    cout << "#### " << i << "th eigen ####" << endl;
    cout << "w=" << w[i] <<endl;
  }
}

//=============================================================================
void dsyev_check_right()
{
  cout << "############ check dsyev right ############" << endl;
  
  srand(unsigned(time(NULL)));
  int N(3);
  
  //// make dsymatrix A  ////
  CPPL::dsymatrix A(N);
  for(int i=0; i<A.n; i++){ for(int j=0; j<A.n; j++){
    A(i,j) =double( rand() /(RAND_MAX/10) );
  }}
  
  //// make wr wi vr ////
  vector<double> w;
  vector<CPPL::dcovector> v;
  
  //// make A_original ////
  CPPL::dsymatrix A_original(A);
  
  //// dsyev ////
  A.dsyev(w,v);
  
  //// print ////
  cout << "A_original=\n" << A_original << endl;
  for(int i=0; i<A.n; i++){
    cout << "#### " << i << "th eigen ####" << endl;
    cout << "w=" << w[i] <<endl;
    cout << "v=\n" << v[i] <<endl;
    cout << "[A]*{x} -lambda*{x} = (Should be zeros)\n"
         << A_original*v[i] -w[i]*v[i] << endl;
  }
}

//=============================================================================
void dsyev_check_left()
{
  cout << "############ check dsyev left ############" << endl;
  
  srand(unsigned(time(NULL)));
  int N(3);
  
  //// make dsymatrix A  ////
  CPPL::dsymatrix A(N);
  for(int i=0; i<A.n; i++){ for(int j=0; j<A.n; j++){
    A(i,j) =double( rand() /(RAND_MAX/10) );
  }}
  
  //// make wr wi vl ////
  vector<double> w;
  vector<CPPL::drovector> v;
  
  //// make A_original ////
  CPPL::dsymatrix A_original(A);
  
  //// dsyev ////
  A.dsyev(w, v);
  
  //// print ////
  cout << "A_original=\n" << A_original << endl;
  for(int i=0; i<A.n; i++){
    cout << "#### " << i << "th eigen ####" << endl;
    cout << "w = " << w[i] << endl;
    cout << "v = " << v[i] << endl;
    cout << "{x}*[A] -{x}*lambda = (Should be zeros)\n"
         << v[i]*A_original -v[i]*w[i] << endl;
  }
}

//=============================================================================
#include <iostream>
#include <cstdlib>
#include <ctime>
#include "cpplapack.h"
using namespace std;

#include "dsysv_check.hpp"
#include "dsyev_check.hpp"
#include "dsygv_check.hpp"

//=============================================================================
/*! main */
int main(int argc, char** argv)
{
  //////// dsysv ////////
  dsysv_check_vector();
  dsysv_check_matrix();
  
  //////// dsyev ////////
  dsyev_check_value();
  dsyev_check_right();
  dsyev_check_left();
  
  //////// dsygv ////////
  dsygv_check_value();
  dsygv_check_both();
  
  return 0;
}