Example Codes of dgematrix LAPACK Functions

//=============================================================================
/*! dgesv_check */
void dgesv_check_vector()
{
  std::cout << "############ check dgesv vector ############" << std::endl;
  
  srand(unsigned(time(NULL)));
  int M(3);
  
  //// make dgematrix A  ////
  CPPL::dgematrix A(M,M);
  for(int i=0; i<A.m; i++){ for(int j=0; j<A.n; j++){
    A(i,j) =double( rand() /(RAND_MAX/10) );
  }}
  
  //// make dcovector y ////
  CPPL::dcovector y(M);
  for(int i=0; i<y.l; i++){
    y(i) =double( rand() /(RAND_MAX/10) );
  }
  
  //// make A_original and y_original ////
  CPPL::dgematrix A_original(A);
  CPPL::dcovector y_original(y);
  std::cout << "A_original=\n" << A_original << std::endl;
  std::cout << "y_original=\n" << y_original << std::endl;
  
  //// solve Ax=y ////
  A.dgesv(y);
    
  //// print A, y and A_original*y ////
  std::cout << "A=\n" << A << std::endl;
  std::cout << "y=\n" << y << std::endl;
  std::cout << "A_original*y=\n" << A_original*y << std::endl;
}

void dgesv_check_matrix()
{
  std::cout << "############ check dgesv matrix ############" << std::endl;
  
  srand(unsigned(time(NULL)));
  int M(3);
  

  //// make dgematrix A  ////
  CPPL::dgematrix A(M,M);
  for(int i=0; i<A.m; i++){ for(int j=0; j<A.n; j++){
    A(i,j) =double( rand() /(RAND_MAX/10) );
  }}
   
  //// make dgematrix Y  ////
  CPPL::dgematrix Y(M,M);
  for(int i=0; i<Y.m; i++){ for(int j=0; j<Y.n; j++){
    Y(i,j) =double( rand() /(RAND_MAX/10) );
  }}
   
  //// make A_original and Y_original ////
  CPPL::dgematrix A_original(A);
  CPPL::dgematrix Y_original(Y);
  std::cout << "A_original=\n" << A_original << std::endl;
  std::cout << "Y_original=\n" << Y_original << std::endl;
  
  //// solve AY=B ////
  A.dgesv(Y);
    
  //// print A, Y and A_original*Y ////
  std::cout << "A=\n" << A << std::endl;
  std::cout << "Y=\n" << Y << std::endl;
  std::cout << "A_original*Y=\n" << A_original*Y << std::endl;
}

//=============================================================================
/*! dgels_check */
void dgels_check_vector()
{
  std::cout << "############ check dgels vector ############" << std::endl;
  
  srand(unsigned(time(NULL)));
  int M(3), N(2);
  
  //// make dgematrix A  ////
  CPPL::dgematrix A(M,N);
  for(int i=0; i<A.m; i++){ for(int j=0; j<A.n; j++){
    A(i,j) =double( rand() /(RAND_MAX/10) );
  }}
  
  //// make dcovector y  ////
  CPPL::dcovector y(M);
  for(int i=0; i<y.l; i++){
    y(i) =double( rand() /(RAND_MAX/10) );
  }
  
  //// make originals ////
  CPPL::dgematrix A_original(A);
  CPPL::dcovector y_original(y);
  
  //// dgels ////
  A.dgels(y);
  
  //// print ////
  std::cout << "######## dgels(vec) ########" << std::endl;
  std::cout << "A_original=\n" << A_original << std::endl;
  std::cout << "y_original=\n" << y_original << std::endl;
  std::cout << "A=\n" << A << std::endl;
  std::cout << "y=\n" << y << std::endl;
  std::cout << "A_original*y=\n" << A_original*y << std::endl;
  
  //// reset ////
  A =A_original;
  y =y_original;
  double r; //residual square
  
  //// dgels ////
  A.dgels(y,r);
  
  //// print ////
  std::cout << "######## dgels(vec, residual) ########" << std::endl;
  std::cout << "A_original=\n" << A_original << std::endl;
  std::cout << "y_original=\n" << y_original << std::endl;
  std::cout << "A=\n" << A << std::endl;
  std::cout << "y=\n" << y << std::endl;
  std::cout << "residual=" << r << std::endl;
  std::cout << "A_original*y=\n" << A_original*y << std::endl;
}

//==============================================================================
void dgels_check_matrix()
{
  std::cout << "############ check dgels matrix ############" << std::endl;
  
  srand(unsigned(time(NULL)));
  int M(3), N(2);
  
  //// make dgematrix A  ////
  CPPL::dgematrix A(M,N);
  for(int i=0; i<A.m; i++){ for(int j=0; j<A.n; j++){
    A(i,j) =double( rand() /(RAND_MAX/10) );
  }}

  //// make dgematrix Y  ////
  CPPL::dgematrix Y(M,M);
  for(int i=0; i<Y.m; i++){ for(int j=0; j<Y.n; j++){
    Y(i,j) =double( rand() /(RAND_MAX/10) );
  }}
   
  //// make originals ////
  CPPL::dgematrix A_original(A);
  CPPL::dgematrix Y_original(Y);
  
  //// dgels ////
  A.dgels(Y);
  
  //// print ////
  std::cout << "######## dgels(mat) ########" << std::endl;
  std::cout << "A_original=\n" << A_original << std::endl;
  std::cout << "Y_original=\n" << Y_original << std::endl;
  std::cout << "A=\n" << A << std::endl;
  std::cout << "Y=\n" << Y << std::endl;
  std::cout << "A_original*Y=\n" << A_original*Y << std::endl;
  
  //// reset ////
  A =A_original;
  Y =Y_original;
  CPPL::drovector R; //residual square
  
  //// dgels ////
  A.dgels(Y,R);
  
  //// print ////
  std::cout << "######## dgels(mat, residual) ########" << std::endl;
  std::cout << "A_original=\n" << A_original << std::endl;
  std::cout << "Y_original=\n" << Y_original << std::endl;
  std::cout << "A=\n" << A << std::endl;
  std::cout << "Y=\n" << Y << std::endl;
  std::cout << "residual=" << R << std::endl;
  std::cout << "A_original*Y=\n" << A_original*Y << std::endl;  
}

//=============================================================================
/*! dgelss_check_vector */
void dgelss_check_vector()
{
  std::cout << "############ check dgelss vector ############" << std::endl;
  
  srand(unsigned(time(NULL)));
  int M(3), N(4);
  int RANK(0);
  double RCOND(-1.0);
  
  //// make dgematrix A  ////
  CPPL::dgematrix A(M,N);
  for(int i=0; i<A.m; i++){ for(int j=0; j<A.n; j++){
    A(i,j) =double( rand() /(RAND_MAX/10) );
  }}

  //// make dcovector b  ////
  CPPL::dcovector b(M);
  for(int i=0; i<b.l; i++){
    b(i) =double( rand() /(RAND_MAX/10) );
  }

  //// make dcovector s  ////
  CPPL::dcovector s;
   
  //// make A_original ////
  CPPL::dgematrix A_original(A);
  
  //// make A_original ////
  CPPL::dcovector b_original(b);
  
  //// dgels ////
  A.dgelss(b,s,RANK,RCOND);
  
  //// print ////
  std::cout << "A_original=\n" << A_original << std::endl;
  std::cout << "b_original=\n" << b_original << std::endl;
  std::cout << "A=\n" << A << std::endl;
  std::cout << "b=\n" << b << std::endl;
  std::cout << "s=\n" << s << std::endl;
  std::cout << "A_original*b=\n" << A_original*b << std::endl;
}


//=============================================================================
/*! dgelss_check_matrix */
void dgelss_check_matrix()
{
  std::cout << "############ check dgelss matrix ############" << std::endl;
  
  srand(unsigned(time(NULL)));
  int M(3), N(4);
  int RANK(0);
  double RCOND(-1.0);
  
  //// make dgematrix A  ////
  CPPL::dgematrix A(M,N);
  for(int i=0; i<A.m; i++){ for(int j=0; j<A.n; j++){
    A(i,j) =double( rand() /(RAND_MAX/10) );
  }}

  //// make dgematrix B  ////
  CPPL::dgematrix b(M,2);
  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 dcovector s  ////
  CPPL::dcovector s;
   
  //// make A_original ////
  CPPL::dgematrix A_original(A);
  
  //// make b_original ////
  CPPL::dgematrix b_original(b);
  
  //// dgels ////
  A.dgelss(b,s,RANK,RCOND);
  
  //// print ////
  std::cout << "A_original=\n" << A_original << std::endl;
  std::cout << "b_original=\n" << b_original << std::endl;
  std::cout << "A=\n" << A << std::endl;
  std::cout << "b=\n" << b << std::endl;
  std::cout << "s=\n" << s << std::endl;
  std::cout << "A_original*b=\n" << A_original*b << std::endl;
}

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

//=============================================================================
void dgeev_check_right()
{
  std::cout << "############ check dgeev right ############" << std::endl;
  
  srand(unsigned(time(NULL)));
  int M(3);
  
  //// make dgematrix A  ////
  CPPL::dgematrix A(M,M);
  for(int i=0; i<A.m; i++){ for(int j=0; j<A.n; j++){
      A(i,j) =double( rand() /(RAND_MAX/10) );
    }}
  
  //// make wr wi vr ////
  std::vector<double> wr, wi;
  std::vector<CPPL::dcovector> vrr, vri;
  
  //// make A_original ////
  CPPL::dgematrix A_original(A);
  
  //// dgeev ////
  A.dgeev(wr, wi ,vrr, vri);
  
  //// print ////
  std::cout << "A_original=\n" << A_original << std::endl;
  for(int i=0; i<A.m; i++){
    std::cout << "#### " << i << "th eigen ####" << std::endl;
    std::cout << "wr=" << wr[i] <<std::endl;
    std::cout << "wi=" << wi[i] <<std::endl;
    std::cout << "vrr=\n" << vrr[i] <<std::endl;
    std::cout << "vri=\n" << vri[i] <<std::endl;
    std::cout << "Real[ [A]*{x} -lambda*{x} ] = (Should be zeros)\n"
              << A_original*vrr[i] -(wr[i]*vrr[i] - wi[i]*vri[i]) 
              << std::endl;
    std::cout << "Imag[ [A]*{x} -lambda*{x} ] = (Should be zeros)\n"
              << A_original*vri[i] -(wr[i]*vri[i] + wi[i]*vrr[i]) 
              << std::endl;
  }
}

//=============================================================================
void dgeev_check_left()
{
  std::cout << "############ check dgeev left ############" << std::endl;
  
  srand(unsigned(time(NULL)));
  int M(3);
  
  //// make dgematrix A  ////
  CPPL::dgematrix A(M,M);
  for(int i=0; i<A.m; i++){ for(int j=0; j<A.n; j++){
      A(i,j) =double( rand() /(RAND_MAX/10) );
    }}
  
  //// make wr wi vl ////
  std::vector<double> wr, wi;
  std::vector<CPPL::drovector> vlr, vli;
  
  //// make A_original ////
  CPPL::dgematrix A_original(A);
  
  //// dgeev ////
  A.dgeev(wr, wi ,vlr, vli);
  
  //// print ////
  std::cout << "A_original=\n" << A_original << std::endl;
  for(int i=0; i<A.m; i++){
    std::cout << "#### " << i << "th eigen ####" << std::endl;
    std::cout << "wr = " << wr[i] << std::endl;
    std::cout << "wi = " << wi[i] << std::endl;
    std::cout << "vlr = " << vlr[i];
    std::cout << "vli = " << vli[i] << std::endl;
    std::cout << "Real[ {x}*[A] -{x}*lambda ] = (Should be zeros)\n"
              << vlr[i]*A_original -(vlr[i]*wr[i] - vli[i]*wi[i]) 
              << std::endl;
    std::cout << "Imag[ {x}*[A] -{x}*lambda ] = (Should be zeros)\n"
              << vli[i]*A_original -(vli[i]*wr[i] + vlr[i]*wi[i]) 
              << std::endl;
  }
}

//=============================================================================
void dggev_check_value()
{
  std::cout << "############ check dggev value ############" << std::endl;
  
  srand(unsigned(time(NULL)));
  int M(3);
  
  //// make dgematrix A and B ////
  CPPL::dgematrix A(M,M), B(M,M);
  for(int i=0; i<A.m; i++){ for(int j=0; j<A.n; j++){
    A(i,j) =double( rand() /(RAND_MAX/10) );
    B(i,j) =double( rand() /(RAND_MAX/10) );
  }}
  
  //// make wr wi vr ////
  std::vector<double> wr, wi;
  
  //// make A_original and B_original ////
  CPPL::dgematrix A_original(A);
  CPPL::dgematrix B_original(B);
  
  //// dggev ////
  A.dggev(B, wr, wi);
  
  //// print ////
  std::cout << "A_original=\n" << A_original << std::endl;
  std::cout << "B_original=\n" << B_original << std::endl;
  for(int i=0; i<A.m; i++){
    std::cout << "#### " << i << "th eigen ####" << std::endl;
    std::cout << "wr=" << wr[i] <<std::endl;
    std::cout << "wi=" << wi[i] <<std::endl;
  }
}

//=============================================================================
void dggev_check_right()
{
  std::cout << "############ check dggev right ############" << std::endl;
  
  srand(unsigned(time(NULL)));
  int M(3);
  
  //// make dgematrix A and B ////
  CPPL::dgematrix A(M,M), B(M,M);
  for(int i=0; i<A.m; i++){ for(int j=0; j<A.n; j++){
    A(i,j) =double( rand() /(RAND_MAX/10) );
    B(i,j) =double( rand() /(RAND_MAX/10) );
  }}
  
  //// make wr wi vr ////
  std::vector<double> wr, wi;
  std::vector<CPPL::dcovector> vrr, vri;
  
  //// make A_original and B_original ////
  CPPL::dgematrix A_original(A);
  CPPL::dgematrix B_original(B);
  
  //// dggev ////
  A.dggev(B, wr, wi ,vrr, vri);
  
  //// print ////
  std::cout << "A_original=\n" << A_original << std::endl;
  std::cout << "B_original=\n" << B_original << std::endl;
  for(int i=0; i<A.m; i++){
    std::cout << "#### " << i << "th eigen ####" << std::endl;
    std::cout << "wr=" << wr[i] <<std::endl;
    std::cout << "wi=" << wi[i] <<std::endl;
    std::cout << "vrr=\n" << vrr[i] <<std::endl;
    std::cout << "vri=\n" << vri[i] <<std::endl;
    std::cout << "Real[ [A]*{x} -lambda*[B]*{x} ] = (Should be zeros)\n"
         << A_original*vrr[i]
      -(wr[i]*B_original*vrr[i] - wi[i]*B_original*vri[i])
         << std::endl;
    std::cout << "Imag[ [A]*{x} -lambda*[B]*{x} ] = (Should be zeros)\n"
         << A_original*vri[i]
      -(wr[i]*B_original*vri[i] + wi[i]*B_original*vrr[i])
         << std::endl;
  }
}

//=============================================================================
void dggev_check_left()
{
  std::cout << "############ check dggev left ############" << std::endl;
  
  srand(unsigned(time(NULL)));
  int M(3);
  
  //// make dgematrix A and B ////
  CPPL::dgematrix A(M,M), B(M,M);
  for(int i=0; i<A.m; i++){ for(int j=0; j<A.n; j++){
    A(i,j) =double( rand() /(RAND_MAX/10) );
    B(i,j) =double( rand() /(RAND_MAX/10) );
  }}
  
  //// make wr wi vl ////
  std::vector<double> wr, wi;
  std::vector<CPPL::drovector> vlr, vli;
  
  //// make A_original and B_original ////
  CPPL::dgematrix A_original(A);
  CPPL::dgematrix B_original(B);
  
  //// dggev ////
  A.dggev(B, wr, wi ,vlr, vli);
  
  //// print ////
  std::cout << "A_original=\n" << A_original << std::endl;
  std::cout << "B_original=\n" << B_original << std::endl;
  for(int i=0; i<A.m; i++){
    std::cout << "#### " << i << "th eigen ####" << std::endl;
    std::cout << "wr = " << wr[i] << std::endl;
    std::cout << "wi = " << wi[i] << std::endl;
    std::cout << "vlr = " << vlr[i];
    std::cout << "vli = " << vli[i] << std::endl;
    std::cout << "Real[ {x}*[A] -lambda*{x}*[B] ] = (Should be zeros)\n"
         << vlr[i]*A_original
      -(wr[i]*vlr[i]*B_original - wi[i]*vli[i]*B_original)
         << std::endl;
    std::cout << "Imag[ [A]*{x} -lambda*{x}*[B] ] = (Should be zeros)\n"
         << vli[i]*A_original
      -(wr[i]*vli[i]*B_original + wi[i]*vlr[i]*B_original)
         << std::endl;
  }
}

//=============================================================================
/*! dgesvd_check */
void dgesvd_check()
{
  std::cout << "############ check dgesvd ############" << std::endl;
  
  srand(unsigned(time(NULL)));
  int M(4), N(3);
  
  //// make dgematrix A  ////
  CPPL::dgematrix A(M,N);
  for(int i=0; i<A.m; i++){ for(int j=0; j<A.n; j++){
  A(i,j) =double( rand() /(RAND_MAX/10) );
  }}
  
  //// make A_original ////
  CPPL::dgematrix A_original(A);
  
  //// make S, U and VT ////
  CPPL::dcovector S;
  CPPL::dgematrix U, VT;
  
  //// SVD A ////
  A.dgesvd(S,U,VT);
  
  //// print A, S, U, and VT ////
  std::cout << "A=\n" << A << std::endl;
  std::cout << "S=\n" << S << std::endl;
  std::cout << "U=\n" << U << std::endl;
  std::cout << "VT=\n" << VT << std::endl;
  
  CPPL::dgematrix S_matrix(M,N);
  S_matrix.zero();
  for(int i=0; i<std::min(M,N); i++){ S_matrix(i,i) =S(i); }
  std::cout << "S_matrix=\n" << S_matrix << std::endl;
  std::cout << "A_original=\n" << A_original << std::endl;
  std::cout << "U*S_matrix*VT=\n" << U*S_matrix*VT << std::endl;
}

//=============================================================================
#include <iostream>
#include <cstdlib>
#include <ctime>
#include "cpplapack.h"

#include "dgesv_check.hpp"
#include "dgels_check.hpp"
#include "dgelss_check.hpp"
#include "dgeev_check.hpp"
#include "dggev_check.hpp"
#include "dgesvd_check.hpp"
#include "dgglse_check.hpp"

//=============================================================================
/*! main */
int main(int argc, char** argv)
{
  //////// dgesv ////////
  dgesv_check_vector();
  dgesv_check_matrix();
  
  //////// dgels ////////
  dgels_check_vector();
  dgels_check_matrix();
  
  //////// dgelss ////////
  dgelss_check_vector();
  dgelss_check_matrix();
  
  //////// dgeev ////////
  dgeev_check_value();
  dgeev_check_right();
  dgeev_check_left();
  
  //////// dggev ////////
  dggev_check_value();
  dggev_check_right();
  dggev_check_left();
  
  //////// dgesvd ////////
  dgesvd_check();

  //////// dgglse ////////
  dgglse_check();
  
  return 0;
}