I am writing this in C and here is a sample section of code. I've compiled this with the same results as my larger section to ensure continuity. More error is seen when the matrix has decimals of higher-precision --I just typed in a simple one to prove my point. I assume this 'error', is in fact not error, and only do to the eigenvector normalization. Any suggestions?
Thanks all!
-------------------------------------------
- Code: Select all
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
void mk_diag(double* diag, double* mat, int n, const int m){
int i;
for(i =0; i < m; ++i)
mat[i*n+i] = diag[i];
}
int main()
{
// --------------------------- test values
double A[] =
{0,1,2,3,
1,2,3,4,
1,3,4,5,
1,1,1,1 };
int n = 4;
// ---------------------------------------------
char jobv_ = 'V';
double *wi,*wr,*U,*work,work_size,*Ui,*D;
int lwork,info;
wi = (double*) malloc( n*sizeof(double));
wr = (double*) malloc( n*sizeof(double));
U = (double*) calloc( n*n,sizeof(double));
D = (double*) calloc( n*n,sizeof(double));
Ui = (double*) calloc( n*n,sizeof(double));
lwork = -1;
dgeev_(&jobv_,&jobv_,&n,A,&n,wr,wi,U,&n,Ui,&n,&work_size,&lwork,&info);
if( info == 0 ) {
lwork = (int)work_size;
work = (double*) calloc( lwork , sizeof( double) );
dgeev_(&jobv_,&jobv_,&n,A,&n,wr,wi,U,&n,Ui,&n,work,&lwork,&info);
if( info == 0 ) {
mk_diag(wr,D, n,n); //create diagonal matrix
} else { return 0; }
free( work );
} else { return 0; }
double *UD, *UDUi, a = 1,b=0;
char ntran_ = 'N';
char ytran_ = 'T';
char ctran_ = 'C';
UD = (double*) malloc( n*n*sizeof(double));
dgemm_( &ntran_,&ntran_,&n,&n,&n,&a,U,&n,D,&n,&b,UD,&n );
UDUi = (double*) malloc( n*n*sizeof(double));
dgemm_( &ntran_,&ytran_,&n,&n,&n,&a,UD,&n,Ui,&n,&b,UDUi,&n );
int i,j;
printf ("\n(A):\n\t");
for (i=0; i<n; ++i) {
for (j=0; j<n; ++j){
printf("[%6.1f] ", UDUi[i*n+j]);
}
printf("\n\t");
}
return 0;
}
