Skip to content
Snippets Groups Projects
Commit 698775c7 authored by joeri.brackenhoff's avatar joeri.brackenhoff
Browse files

3D

parents 0efc1592 8c64dbf4
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 955 additions and 20 deletions
.gitignore
+ 17
− 0
View file @ 698775c7
*.[oa] *.[oa]
*.eps
bin/* bin/*
*.su *.su
*.bin *.bin
...@@ -20,6 +21,7 @@ utils/makemod ...@@ -20,6 +21,7 @@ utils/makemod
utils/makewave utils/makewave
utils/mat2su utils/mat2su
utils/syn2d utils/syn2d
utils/green3D
fdelmodc/fdelmodc fdelmodc/fdelmodc
include/genfft.h include/genfft.h
Make_include Make_include
...@@ -29,3 +31,18 @@ marchenko_full/fmute ...@@ -29,3 +31,18 @@ marchenko_full/fmute
marchenko_full/marchenko marchenko_full/marchenko
corrvir/corrvir corrvir/corrvir
fdemmodc/fdemmodc fdemmodc/fdemmodc
FFTlib/test/cc1test
FFTlib/test/cc2dtest
FFTlib/test/ccmtest
FFTlib/test/rc1test
FFTlib/test/rc2dtest
FFTlib/test/rcmtest
marchenko_applications/HomG
marchenko_applications/MuteSnap
marchenko_applications/combine
marchenko_applications/gmshift
marchenko_applications/iba
marchenko_applications/marchenko_app
marchenko_applications/reshape_su
raytime3d/raytime3d
MDD/mdd
FFTlib/cc1fft.c
+ 30
− 0
View file @ 698775c7
#include "genfft.h" #include "genfft.h"
#ifdef MKL
#include "mkl_dfti.h"
void dfti_status_print(MKL_LONG status);
#endif
/** /**
* NAME: cc1fft * NAME: cc1fft
...@@ -56,6 +60,10 @@ void cc1fft(complex *data, int n, int sign) ...@@ -56,6 +60,10 @@ void cc1fft(complex *data, int n, int sign)
static complex *work; static complex *work;
REAL scl; REAL scl;
complex *y; complex *y;
#elif defined(MKL)
static DFTI_DESCRIPTOR_HANDLE handle=0;
static int nprev=0;
MKL_LONG Status;
#endif #endif
#if defined(HAVE_LIBSCS) #if defined(HAVE_LIBSCS)
...@@ -90,6 +98,28 @@ void cc1fft(complex *data, int n, int sign) ...@@ -90,6 +98,28 @@ void cc1fft(complex *data, int n, int sign)
nprev = n; nprev = n;
} }
acmlcc1fft(sign, scl, inpl, n, data, 1, y, 1, work, &isys); acmlcc1fft(sign, scl, inpl, n, data, 1, y, 1, work, &isys);
#elif defined(MKL)
if (n != nprev) {
DftiFreeDescriptor(&handle);
Status = DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_COMPLEX, 1, (MKL_LONG)n);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCreateDescriptor FAIL\n");
}
Status = DftiCommitDescriptor(handle);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCommitDescriptor FAIL\n");
}
nprev = n;
}
if (sign < 0) {
Status = DftiComputeBackward(handle, data);
}
else {
Status = DftiComputeForward(handle, data);
}
#else #else
cc1_fft(data, n, sign); cc1_fft(data, n, sign);
#endif #endif
......
FFTlib/cc2dfft.c
+ 32
− 0
View file @ 698775c7
#include "genfft.h" #include "genfft.h"
#ifdef MKL
#include "mkl_dfti.h"
void dfti_status_print(MKL_LONG status);
#endif
/** /**
* NAME: cc2dfft * NAME: cc2dfft
...@@ -60,6 +64,10 @@ void cc2dfft(complex *data, int nx, int ny, int ldx, int sign) ...@@ -60,6 +64,10 @@ void cc2dfft(complex *data, int nx, int ny, int ldx, int sign)
REAL scl; REAL scl;
complex *y; complex *y;
complex *tmp; complex *tmp;
#elif defined(MKL)
static DFTI_DESCRIPTOR_HANDLE handle=0;
static int nxprev=0, nyprev=0;
MKL_LONG Status, N[2];
#else #else
int i, j; int i, j;
complex *tmp; complex *tmp;
...@@ -110,6 +118,30 @@ void cc2dfft(complex *data, int nx, int ny, int ldx, int sign) ...@@ -110,6 +118,30 @@ void cc2dfft(complex *data, int nx, int ny, int ldx, int sign)
else { else {
cc1fft(data, nx, sign); cc1fft(data, nx, sign);
} }
#elif defined(MKL)
if (nx != nxprev || ny != nyprev) {
DftiFreeDescriptor(&handle);
N[0] = nx; N[1] = ny;
Status = DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_COMPLEX, 2, N);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCreateDescriptor FAIL\n");
}
Status = DftiCommitDescriptor(handle);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCommitDescriptor FAIL\n");
}
nxprev = nx;
nyprev = ny;
}
if (sign < 0) {
Status = DftiComputeBackward(handle, data);
}
else {
Status = DftiComputeForward(handle, data);
}
#else #else
if (ny != 1) { if (ny != 1) {
ccmfft(data, nx, ny, ldx, sign); ccmfft(data, nx, ny, ldx, sign);
......
FFTlib/ccmfft.c
+ 35
− 0
View file @ 698775c7
#include "genfft.h" #include "genfft.h"
#ifdef MKL
#include "mkl_dfti.h"
void dfti_status_print(MKL_LONG status);
#endif
/** /**
* NAME: ccmfft * NAME: ccmfft
...@@ -59,6 +63,11 @@ void ccmfft(complex *data, int n1, int n2, int ld1, int sign) ...@@ -59,6 +63,11 @@ void ccmfft(complex *data, int n1, int n2, int ld1, int sign)
static complex *work; static complex *work;
REAL scl; REAL scl;
complex *y; complex *y;
#elif defined(MKL)
static DFTI_DESCRIPTOR_HANDLE handle=0;
static int nprev=0;
MKL_LONG Status;
int j;
#endif #endif
#if defined(HAVE_LIBSCS) #if defined(HAVE_LIBSCS)
...@@ -89,6 +98,32 @@ void ccmfft(complex *data, int n1, int n2, int ld1, int sign) ...@@ -89,6 +98,32 @@ void ccmfft(complex *data, int n1, int n2, int ld1, int sign)
nprev = n1; nprev = n1;
} }
acmlccmfft(sign, scl, inpl, n2, n1, data, 1, ld1, y, 1, ld1, work, &isys); acmlccmfft(sign, scl, inpl, n2, n1, data, 1, ld1, y, 1, ld1, work, &isys);
#elif defined(MKL)
if (n1 != nprev) {
DftiFreeDescriptor(&handle);
Status = DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_COMPLEX, 1, (MKL_LONG)n1);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCreateDescriptor FAIL\n");
}
Status = DftiCommitDescriptor(handle);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCommitDescriptor FAIL\n");
}
nprev = n1;
}
if (sign < 0) {
for (j=0; j<n2; j++) {
Status = DftiComputeBackward(handle, &data[j*ld1]);
}
}
else {
for (j=0; j<n2; j++) {
Status = DftiComputeForward(handle, &data[j*ld1]);
}
}
#else #else
ccm_fft(data, n1, n2, ld1, sign); ccm_fft(data, n1, n2, ld1, sign);
#endif #endif
......
FFTlib/cr1fft.c
+ 51
− 4
View file @ 698775c7
#include "genfft.h" #include "genfft.h"
#ifdef MKL
#include "mkl_dfti.h"
void dfti_status_print(MKL_LONG status);
#endif
/** /**
* NAME: cr1fft * NAME: cr1fft
...@@ -52,10 +56,12 @@ void cr1fft(complex *cdata, REAL *rdata, int n, int sign) ...@@ -52,10 +56,12 @@ void cr1fft(complex *cdata, REAL *rdata, int n, int sign)
static int isys; static int isys;
static REAL *work, *table, scale=1.0; static REAL *work, *table, scale=1.0;
REAL scl; REAL scl;
#elif defined(FFTW3) #elif defined(MKL)
static int nprev=0; static DFTI_DESCRIPTOR_HANDLE handle=0;
int iopt, ier; static int nprev=0;
static float *work; REAL *tmp;
MKL_LONG Status;
int i;
#endif #endif
#if defined(HAVE_LIBSCS) #if defined(HAVE_LIBSCS)
...@@ -91,6 +97,47 @@ void cr1fft(complex *cdata, REAL *rdata, int n, int sign) ...@@ -91,6 +97,47 @@ void cr1fft(complex *cdata, REAL *rdata, int n, int sign)
nprev = n; nprev = n;
} }
acmlcrfft(one, n, rdata, work, &isys); acmlcrfft(one, n, rdata, work, &isys);
#elif defined(MKL)
if (n != nprev) {
DftiFreeDescriptor(&handle);
Status = DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_REAL, 1, (MKL_LONG)n);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCreateDescriptor FAIL\n");
}
Status = DftiSetValue(handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
Status = DftiSetValue(handle, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
if (! DftiErrorClass(Status, DFTI_NO_ERROR)) {
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
Status = DftiCommitDescriptor(handle);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCommitDescriptor FAIL\n");
}
nprev = n;
}
tmp = (float *)malloc(n*sizeof(float));
Status = DftiComputeBackward(handle, (MKL_Complex8 *)cdata, tmp);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiComputeBackward FAIL\n");
}
rdata[0] = tmp[0];
if (sign < 0) {
for (i=1; i<n; i++) rdata[i] = -sign*tmp[n-i];
}
else {
for (i=1; i<n; i++) rdata[i] = tmp[i];
}
free(tmp);
#else #else
cr1_fft(cdata, rdata, n, sign); cr1_fft(cdata, rdata, n, sign);
#endif #endif
......
FFTlib/crmfft.c
+ 55
− 0
View file @ 698775c7
#include "genfft.h" #include "genfft.h"
#include <string.h> #include <string.h>
#ifdef MKL
#include "mkl_dfti.h"
void dfti_status_print(MKL_LONG status);
#endif
/** /**
* NAME: crmfft * NAME: crmfft
...@@ -66,6 +70,12 @@ void crmfft(complex *cdata, REAL *rdata, int n1, int n2, int ldc, int ldr, int s ...@@ -66,6 +70,12 @@ void crmfft(complex *cdata, REAL *rdata, int n1, int n2, int ldc, int ldr, int s
static int isys; static int isys;
static REAL *work; static REAL *work;
REAL scl, *data; REAL scl, *data;
#elif defined(MKL)
static DFTI_DESCRIPTOR_HANDLE handle=0;
static int nprev=0;
REAL *tmp;
MKL_LONG Status;
int i, j;
#endif #endif
#if defined(HAVE_LIBSCS) #if defined(HAVE_LIBSCS)
...@@ -127,6 +137,51 @@ void crmfft(complex *cdata, REAL *rdata, int n1, int n2, int ldc, int ldr, int s ...@@ -127,6 +137,51 @@ void crmfft(complex *cdata, REAL *rdata, int n1, int n2, int ldc, int ldr, int s
for (j=0; j<n2; j++) { for (j=0; j<n2; j++) {
memcpy(&rdata[j*ldr],&data[j*n1],n1*sizeof(REAL)); memcpy(&rdata[j*ldr],&data[j*n1],n1*sizeof(REAL));
} }
#elif defined(MKL)
if (n1 != nprev) {
DftiFreeDescriptor(&handle);
Status = DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_REAL, 1, (MKL_LONG)n1);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCreateDescriptor FAIL\n");
}
Status = DftiSetValue(handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
Status = DftiSetValue(handle, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
//This options is what we would like, but is depreciated in the future
//Status = DftiSetValue(handle, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_REAL);
if (! DftiErrorClass(Status, DFTI_NO_ERROR)) {
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
Status = DftiCommitDescriptor(handle);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCommitDescriptor FAIL\n");
}
nprev = n1;
}
tmp = (float *)malloc(n1*sizeof(float));
for (j=0; j<n2; j++) {
Status = DftiComputeBackward(handle, (MKL_Complex8 *)&cdata[j*ldc], tmp);
rdata[j*ldr] = tmp[0];
if (sign < 0) {
for (i=1; i<n1; i++) rdata[j*ldr+i] = -sign*tmp[n1-i];
}
else {
for (i=1; i<n1; i++) rdata[j*ldr+i] = tmp[i];
}
}
free(tmp);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiComputeBackward FAIL\n");
}
#else #else
crm_fft(cdata, rdata, n1, n2, ldc, ldr, sign); crm_fft(cdata, rdata, n1, n2, ldc, ldr, sign);
#endif #endif
......
FFTlib/rc1fft.c
+ 64
− 0
View file @ 698775c7
#include "genfft.h" #include "genfft.h"
#include <string.h> #include <string.h>
#ifdef MKL
#include "mkl_dfti.h"
void dfti_status_print(MKL_LONG status);
#endif
/** /**
* NAME: rc1fft * NAME: rc1fft
...@@ -52,6 +56,11 @@ void rc1fft(REAL *rdata, complex *cdata, int n, int sign) ...@@ -52,6 +56,11 @@ void rc1fft(REAL *rdata, complex *cdata, int n, int sign)
static int isys; static int isys;
static REAL *work, *table, scale=1.0; static REAL *work, *table, scale=1.0;
REAL scl, *data; REAL scl, *data;
#elif defined(MKL)
static DFTI_DESCRIPTOR_HANDLE handle=0;
static int nprev=0;
MKL_LONG Status;
int i;
#endif #endif
#if defined(HAVE_LIBSCS) #if defined(HAVE_LIBSCS)
...@@ -92,6 +101,49 @@ void rc1fft(REAL *rdata, complex *cdata, int n, int sign) ...@@ -92,6 +101,49 @@ void rc1fft(REAL *rdata, complex *cdata, int n, int sign)
} }
cdata[n/2].i=0.0; cdata[n/2].i=0.0;
free(data); free(data);
#elif defined(MKL)
if (n != nprev) {
DftiFreeDescriptor(&handle);
Status = DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_REAL, 1, (MKL_LONG)n);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCreateDescriptor FAIL\n");
}
Status = DftiSetValue(handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
/*
Status = DftiSetValue(handle, DFTI_FORWARD_DOMAIN, DFTI_REAL);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
*/
Status = DftiSetValue(handle, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
if (! DftiErrorClass(Status, DFTI_NO_ERROR)) {
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
Status = DftiCommitDescriptor(handle);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCommitDescriptor FAIL\n");
}
nprev = n;
}
Status = DftiComputeForward(handle, rdata, (MKL_Complex8 *)cdata);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiComputeForward FAIL\n");
}
for (i=1; i<((n-1)/2)+1; i++) {
cdata[i].i *= -sign;
}
#else #else
rc1_fft(rdata, cdata, n, sign); rc1_fft(rdata, cdata, n, sign);
#endif #endif
...@@ -99,6 +151,18 @@ void rc1fft(REAL *rdata, complex *cdata, int n, int sign) ...@@ -99,6 +151,18 @@ void rc1fft(REAL *rdata, complex *cdata, int n, int sign)
return; return;
} }
#ifdef MKL
void dfti_status_print(MKL_LONG status)
{
MKL_LONG class_error;
char* error_message;
error_message = DftiErrorMessage(status);
printf("error_message = %s \n", error_message);
return;
}
#endif
/****************** NO COMPLEX DEFINED ******************/ /****************** NO COMPLEX DEFINED ******************/
......
FFTlib/rcmfft.c
+ 47
− 0
View file @ 698775c7
#include "genfft.h" #include "genfft.h"
#include <string.h> #include <string.h>
#ifdef MKL
#include "mkl_dfti.h"
void dfti_status_print(MKL_LONG status);
#endif
/** /**
* NAME: rcmfft * NAME: rcmfft
...@@ -64,6 +68,11 @@ void rcmfft(REAL *rdata, complex *cdata, int n1, int n2, int ldr, int ldc, int s ...@@ -64,6 +68,11 @@ void rcmfft(REAL *rdata, complex *cdata, int n1, int n2, int ldr, int ldc, int s
static int isys; static int isys;
static REAL *work; static REAL *work;
REAL scl, *data; REAL scl, *data;
#elif defined(MKL)
static DFTI_DESCRIPTOR_HANDLE handle=0;
static int nprev=0;
MKL_LONG Status;
int i,j;
#endif #endif
#if defined(HAVE_LIBSCS) #if defined(HAVE_LIBSCS)
...@@ -120,6 +129,44 @@ void rcmfft(REAL *rdata, complex *cdata, int n1, int n2, int ldr, int ldc, int s ...@@ -120,6 +129,44 @@ void rcmfft(REAL *rdata, complex *cdata, int n1, int n2, int ldr, int ldc, int s
cdata[j*ldc+n1/2].i=0.0; cdata[j*ldc+n1/2].i=0.0;
} }
free(data); free(data);
#elif defined(MKL)
if (n1 != nprev) {
DftiFreeDescriptor(&handle);
Status = DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_REAL, 1, (MKL_LONG)n1);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCreateDescriptor FAIL\n");
}
Status = DftiSetValue(handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
Status = DftiSetValue(handle, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
if (! DftiErrorClass(Status, DFTI_NO_ERROR)) {
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
Status = DftiCommitDescriptor(handle);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCommitDescriptor FAIL\n");
}
nprev = n1;
}
Status = DftiComputeForward(handle, rdata, (MKL_Complex8 *)cdata);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiComputeForward FAIL\n");
}
for (j=0; j<n2; j++) {
Status = DftiComputeForward(handle, &rdata[j*ldr], (MKL_Complex8 *)&cdata[j*ldc]);
for (i=1; i<((n1-1)/2)+1; i++) {
cdata[j*ldc+i].i *= -sign;
}
}
#else #else
rcm_fft(rdata, cdata, n1, n2, ldr, ldc, sign); rcm_fft(rdata, cdata, n1, n2, ldr, ldc, sign);
#endif #endif
......
FFTlib/test/Makefile
+ 1
− 1
View file @ 698775c7
...@@ -2,7 +2,7 @@ ...@@ -2,7 +2,7 @@
include ../../Make_include include ../../Make_include
LIBS += -L$L -lgenfft -lm LIBS += -L$L -lgenfft $(FFT) -lm
ALL: cc1test cc2dtest ccmtest rc1test rc2dtest rcmtest rc1loop cc1loop ALL: cc1test cc2dtest ccmtest rc1test rc2dtest rcmtest rc1loop cc1loop
......
FFTlib/test/cc1loop.c
+ 2
− 2
View file @ 698775c7
#include <genfft.h> #include <genfft.h>
#include <time.h> #include <time.h>
#include <kiss_fft.h>
main () { int main () {
int j,i,n,sign, isign; int j,i,n,sign, isign;
int N, Nmax=600, Nitcc; int N, Nmax=600, Nitcc;
...@@ -65,5 +64,6 @@ main () { ...@@ -65,5 +64,6 @@ main () {
N += 1; N += 1;
} }
return 0;
} }
FFTlib/test/cc1test.c
+ 2
− 1
View file @ 698775c7
#include <genfft.h> #include <genfft.h>
#include <time.h> #include <time.h>
main () { int main () {
int j,i,n,sign, isign; int j,i,n,sign, isign;
int N, Nmax=8192, Nitcc; int N, Nmax=8192, Nitcc;
...@@ -79,5 +79,6 @@ main () { ...@@ -79,5 +79,6 @@ main () {
k += 1.0; k += 1.0;
} }
return 0;
} }
FFTlib/test/cc2dtest.c
+ 2
− 1
View file @ 698775c7
#include <genfft.h> #include <genfft.h>
#include <time.h> #include <time.h>
main () { int main () {
int j,i,n,sign, isign; int j,i,n,sign, isign;
int N, Nmax=1024, Nitcc, Nlot=1024, ld1; int N, Nmax=1024, Nitcc, Nlot=1024, ld1;
...@@ -142,5 +142,6 @@ main () { ...@@ -142,5 +142,6 @@ main () {
Nlot *= 2; Nlot *= 2;
} }
return 0;
} }
FFTlib/test/ccmtest.c
+ 2
− 1
View file @ 698775c7
#include <genfft.h> #include <genfft.h>
#include <time.h> #include <time.h>
main () { int main () {
int j,i,n,sign, isign; int j,i,n,sign, isign;
int N, Nmax=4097, Nitcc, Nlot=513, ld1; int N, Nmax=4097, Nitcc, Nlot=513, ld1;
...@@ -128,5 +128,6 @@ main () { ...@@ -128,5 +128,6 @@ main () {
k += 1.0; k += 1.0;
} }
return 0;
} }
FFTlib/test/rc1loop.c
+ 2
− 1
View file @ 698775c7
...@@ -5,7 +5,7 @@ void crdft(complex *cdata, REAL *rdata, int n, int sign); ...@@ -5,7 +5,7 @@ void crdft(complex *cdata, REAL *rdata, int n, int sign);
void rcdft(REAL *rdata, complex *cdata, int n, int sign); void rcdft(REAL *rdata, complex *cdata, int n, int sign);
main () { int main () {
int j,i,n,k,sign, isign; int j,i,n,k,sign, isign;
int N, Nmax=600, Nitcc; int N, Nmax=600, Nitcc;
...@@ -94,5 +94,6 @@ main () { ...@@ -94,5 +94,6 @@ main () {
N += 1; N += 1;
} }
return 0;
} }
FFTlib/test/rc1test.c
+ 23
− 6
View file @ 698775c7
#include <genfft.h> #include <genfft.h>
#include <time.h> #include <time.h>
main () { int main () {
int j,i,n,sign, isign; int j,i,n,sign, isign;
int N, Nmax=8192, Nitcc; int N, Nmax=8192, Nitcc;
...@@ -24,25 +24,41 @@ main () { ...@@ -24,25 +24,41 @@ main () {
N = 16; N = 16;
k = 5.0; k = 5.0;
sign = 1; sign = -1;
isign = -1; isign = 1;
while (N <= Nmax) { while (N <= Nmax) {
/* Initialize the data */ /* Initialize the data */
for (i=0;i<N;i++) { for (i=0;i<N;i++) {
c_data[i] = (float)-0.1+0.5*(N/2-i); c_data[i] = (float)-0.1+0.5*(N/3-i)*sin(i*M_PI/N);
// c_data[i] = 0.0; // c_data[i] = 0.0;
} }
// c_data[0] = 1.0; // c_data[0] = 1.0;
t = 0.0; t = 0.0;
/*
N=16;
scl = 1.0/(float)N;
for (j=0;j<N;j++) data[j] = c_data[j];
rc1fft(data, cdata, N, sign);
cr1fft(cdata, data, N, isign);
for (i=0; i<N; i++)
fprintf(stderr,"%s: i = %d data = %f Ref-data = %f Complex=%f,%f\n", machine, i, data[i]*scl, c_data[i],cdata[i/2].r,cdata[i/2].i);
for (j=0;j<N;j++) data[j] = c_data[j];
rc1_fft(data, cdata, N, sign);
cr1_fft(cdata, data, N, isign);
for (i=0; i<N; i++)
fprintf(stderr,"%s: i = %d data = %f Ref-data = %f Complex=%f,%f\n", machine, i, data[i]*scl, c_data[i],cdata[i/2].r,cdata[i/2].i);
return;
*/
/* FFT */ /* FFT */
for (i=0; i<2500; i++) { for (i=0; i<2500; i++) {
for (j=0;j<N;j++) data[j] = c_data[j]; for (j=0;j<N;j++) data[j] = c_data[j];
t0 = wallclock_time(); t0 = wallclock_time();
rc1fft(data, cdata, N, sign); rc1fft(data, cdata, N, sign);
cr1fft(cdata, data, N, isign); cr1fft(cdata, data, N, isign);
// cr1_fft(cdata, data, N, isign);
t1 = wallclock_time(); t1 = wallclock_time();
t += t1-t0; t += t1-t0;
} }
...@@ -51,9 +67,9 @@ main () { ...@@ -51,9 +67,9 @@ main () {
scl = 1.0/(float)N; scl = 1.0/(float)N;
for (i=0; i<N; i++) { for (i=0; i<N; i++) {
/* /*
fprintf(stderr,"%s: i = %d data = %f C-data = %f\n", machine, i, data[i].r*scl, c_data[i].r); fprintf(stderr,"%s: i = %d data = %f Ref-data = %f Complex=%f,%f\n", machine, i, data[i]*scl, c_data[i],cdata[i/2].r,cdata[i/2].i);
fprintf(stderr,"%s: i = %d data = %f C-data = %f\n", machine, i, data[i].i*scl, c_data[i].i);
*/ */
if (c_data[i] != 0.0) { if (c_data[i] != 0.0) {
diff = fabs((data[i]*scl - c_data[i]) / c_data[i]); diff = fabs((data[i]*scl - c_data[i]) / c_data[i]);
} }
...@@ -77,5 +93,6 @@ main () { ...@@ -77,5 +93,6 @@ main () {
k += 1.0; k += 1.0;
} }
return 0;
} }
FFTlib/test/rc2dtest.c
+ 2
− 2
View file @ 698775c7
#include <genfft.h> #include <genfft.h>
#include <time.h> #include <time.h>
main () { int main () {
int l,j,i,n,sign, isign, ldr, ldc; int l,j,i,n,sign, isign, ldr, ldc;
int N, Nmax=2048, Nlot=150, Nitcc; int N, Nmax=2048, Nlot=150, Nitcc;
...@@ -141,6 +141,6 @@ main () { ...@@ -141,6 +141,6 @@ main () {
Nlot *= 2; Nlot *= 2;
} }
return 0;
} }
FFTlib/test/rcmtest.c
+ 2
− 1
View file @ 698775c7
#include <genfft.h> #include <genfft.h>
#include <time.h> #include <time.h>
main () { int main () {
int l,j,i,n,sign, isign, ldr, ldc; int l,j,i,n,sign, isign, ldr, ldc;
int N, Nmax=8193, Nlot=150, Nitcc; int N, Nmax=8193, Nlot=150, Nitcc;
...@@ -82,5 +82,6 @@ main () { ...@@ -82,5 +82,6 @@ main () {
k += 1.0; k += 1.0;
} }
return 0;
} }
MDD/Makefile 0 → 100644
+ 61
− 0
View file @ 698775c7
# Makefile
include ../Make_include
########################################################################
# define general include and system library
ALLINC = -I.
#BLAS libs with Intel compiler
#LIBS += -mkl -L$L -lgenfft $(LIBSM)
#General BLAS library
#LIBS += -L$L -lgenfft $(LIBSM)
#General BLAS library
#LIBS += $(BLAS)
#CFLAGS += -I$(MKLROOT)/include
#LIBS += -lblas -llapack -L$L -lgenfft $(LIBSM) -lc -lm
all: mdd
PRG = mdd
SRCC = $(PRG).c \
atopkge.c \
docpkge.c \
getpars.c \
readShotData.c \
writeEigen.c \
deconvolve.c \
computeMatrixInverse.c \
getFileInfo.c \
verbosepkg.c \
name_ext.c \
wallclock_time.c
OBJC = $(SRCC:%.c=%.o)
$(PRG): $(OBJC)
$(CC) $(LDFLAGS) $(CFLAGS) $(OPTC) -o $(PRG) $(OBJC) $(LIBS)
install: $(PRG)
cp $(PRG) $B
clean:
rm -f core $(OBJC) $(OBJM) $(PRG)
realclean:
rm -f core $(OBJC) $(OBJM) $(PRG) $B/$(PRG)
print: Makefile $(SRC)
$(PRINT) $?
@touch print
count:
@wc $(SRC)
tar:
@tar cf $(PRG).tar Makefile $(SRC) && compress $(PRG).tar
MDD/atopkge.c 0 → 120000
+ 1
− 0
View file @ 698775c7
../utils/atopkge.c
\ No newline at end of file
MDD/computeMatrixInverse.c 0 → 100644
+ 524
− 0
View file @ 698775c7
#include<math.h>
#include<stdlib.h>
#include<stdio.h>
#include <assert.h>
#define MAX(x,y) ((x) > (y) ? (x) : (y))
/* Cholesky based inverse */
void cpotrf_(char *uplo, int *N, float *A, int *lda, int *info);
void cpotri_(char *uplo, int *N, float *A, int *lda, int *info);
/* LU based inverse */
void cgetrf_(int *M, int *N, float *A, int *lda, int *ipvt, int *info);
void cgetri_(int *N, float *A, int *lda, int *ipvt, float *work, int *lwork, int *info);
void zgetrf_(int *M, int *N, double *A, int *lda, int *ipvt, int *info);
void zgetri_(int *N, double *A, int *lda, int *ipvt, double *work, int *lwork, int *info);
int ilaenv_(int *ispec, char *name, char *opts, int *n1, int *n2, int *n3, int *n4);
/* SVD based inverse */
void cgesvd_(char *jobu, char *jobvt, int *M, int *N, float *A, int *lda, float *S, float *U, int *ldu, float *vt, int *ldvt, float *work, int *lwork, float *rwork, int *info);
void zgesvd_(char *jobu, char *jobvt, int *M, int *N, double *A, int *lda, double *S, double *U, int *ldu, double *vt, int *ldvt, double *work, int *lwork, double *rwork, int *info);
void cgesdd_(char *jobz, int *M, int *N, float *A, int *lda, float *S, float *U, int *ldu, float *vt, int *ldvt, float *work, int *lwork, float *rwork, int *iwork, int *info);
/* Eigenvalues */
void zgeev_(char *jobvl, char *jobvr, int *N, double *A, int *lda, double *S, double *vl, int *ldvl, double *vr, int *ldvr,
double *work, int *lwork, double *rwork, int *info);
typedef struct { /* complex number */
float r,i;
} complex;
void computeMatrixInverse(complex *matrix, int nxm, int rthm, float eps_a, float eps_r, float numacc, int eigenvalues, float *eigen, int iw, int verbose)
{
int i,j,k,N,lda,info,lwork,*ipvt;
float energy;
complex tmp, one, *work;
char *uplo;
uplo = "U";
lda = N = nxm;
one.r = 1.0;
one.i = 0.0;
if (rthm==0) {
energy=0.0;
if (eps_r != 0.0) {
for (i=0; i<nxm; i++) {
for (j=0; j<nxm; j++) {
tmp = matrix[i*nxm+j];
energy += sqrt(tmp.r*tmp.r+tmp.i*tmp.i);
}
// fprintf(stderr,"i=%d energy=%e\n", i, energy);
}
}
if (verbose>1) fprintf(stderr,"energy=%e eps_r=%e eps_a=%e\n", energy, eps_r*energy, eps_a);
/* add small value at diagonal */
#pragma ivdep
for (i=0; i<nxm; i++) {
tmp.r = eps_r*energy+eps_a;
matrix[i*nxm+i].r+=tmp.r;
}
/* Cholesky based matrix inversion */
cpotrf_(uplo, &N, &matrix[0].r, &lda, &info);
assert (info == 0);
cpotri_(uplo, &N, &matrix[0].r, &lda, &info);
assert (info == 0);
/* fill lower part of inverse matrix */
for (i=0; i<nxm; i++) {
#pragma ivdep
for (j=i+1; j<nxm; j++) {
matrix[i*nxm+j].r=matrix[j*nxm+i].r;
matrix[i*nxm+j].i=-1.0*matrix[j*nxm+i].i;
}
}
}
else if (rthm==1) {
int ispec, n1, nb;
char *name , *opts;
ispec = 1;
name = "CGETRI";
n1 = nxm;
nb = ilaenv_(&ispec, name, opts, &n1, &n1, &n1, &n1);
nb = MAX(1,nb);
lwork = nb*nxm;
ipvt = (int *)malloc(nxm*sizeof(int));
work = (complex *)malloc(lwork*sizeof(complex));
energy=0.0;
if (eps_r != 0.0) {
for (i=0; i<nxm; i++) {
for (j=0; j<nxm; j++) {
tmp = matrix[i*nxm+j];
energy += sqrt(tmp.r*tmp.r+tmp.i*tmp.i);
}
}
}
if (verbose>1) fprintf(stderr,"eps_r=%e eps_a=%e\n", eps_r*energy, eps_a);
/* add small value at diagonal */
for (i=0; i<nxm; i++) {
tmp.r = eps_r*energy+eps_a;
matrix[i*nxm+i].r+=tmp.r;
}
/* LU based matrix inversion */
cgetrf_(&nxm, &nxm, &matrix[0].r, &nxm, ipvt, &info);
assert (info == 0);
cgetri_(&nxm, &matrix[0].r, &nxm, ipvt, &work[0].r, &lwork, &info);
assert (info == 0);
free(ipvt);
free(work);
}
else if (rthm==2) { /* SVD general algorithm most accurate */
float *rwork, *S;
double S0,Si;
complex *U, *VT, a, b;
char *jobu, *jobvt;
int neig;
energy=0.0;
if (eps_r != 0.0) {
for (i=0; i<nxm; i++) {
for (j=0; j<nxm; j++) {
tmp = matrix[i*nxm+j];
energy += sqrt(tmp.r*tmp.r+tmp.i*tmp.i);
}
}
fprintf(stderr,"energy = %e\n", energy);
}
if (verbose>1) fprintf(stderr,"eps_r=%e eps_a=%e\n", eps_r*energy, eps_a);
/* add small value at diagonal */
for (i=0; i<nxm; i++) {
tmp.r = eps_r*energy+eps_a;
matrix[i*nxm+i].r+=tmp.r;
}
jobu = "A";
jobvt = "A";
lda = N = nxm;
lwork = N*8;
S = (float *)malloc(N*sizeof(float));
U = (complex *)malloc(N*N*sizeof(complex));
VT = (complex *)malloc(N*N*sizeof(complex));
work = (complex *)malloc(lwork*sizeof(complex));
rwork = (float *)malloc(5*N*sizeof(float));
/* Compute SVD */
cgesvd_(jobu, jobvt, &N, &N, &matrix[0].r, &lda, S, &U[0].r, &lda, &VT[0].r,
&lda, &work[0].r, &lwork, rwork, &info);
assert (info == 0);
if (eigenvalues) {
for (i=0; i<N; i++) {
eigen[i] = S[i];
}
}
/* Compute inverse */
S0 = S[0];
neig = 0;
for (i=0; i<N; i++) {
/* fprintf(stderr,"S[%d] = %e ",i,S[i]);*/
Si = S[i];
if ((Si/S0) > numacc) { S[i]=1.0/S[i]; neig++; }
else S[i] = 0.0;
/*S[i]=1.0/(S[i]+eps_r*S[0]);*/
/* fprintf(stderr,"S^-1[%d] = %e\n",i,S[i]);*/
}
if(verbose) fprintf(stderr,"fraction of eigenvalues used = %.3f\n",(float)(neig/((float)N)));
for (j=0; j<N; j++) {
for (i=0; i<N; i++) {
U[j*N+i].r=S[j]*U[j*N+i].r;
U[j*N+i].i=-1.0*S[j]*U[j*N+i].i;
}
}
for (j=0; j<N; j++) {
for (i=0; i<N; i++) {
tmp.r = tmp.i = 0.0;
for (k=0; k<N; k++) {
a = U[k*N+j];
b.r = VT[i*N+k].r;
b.i = -1.0*VT[i*N+k].i;
tmp.r += (a.r*b.r-a.i*b.i);
tmp.i += (a.r*b.i+a.i*b.r);
}
matrix[j*nxm+i] = tmp;
}
}
free(U);
free(VT);
free(S);
free(work);
free(rwork);
}
else if (rthm==3) { /* SVD algorithm Divide and Conquerer less accurate */
/* CGESDD*/
int *iwork;
int neig;
float *rwork, *S;
double S0,Si;
complex *U, *VT, a, b;
char *jobz;
energy=0.0;
if (eps_r != 0.0) {
for (i=0; i<nxm; i++) {
for (j=0; j<nxm; j++) {
tmp = matrix[i*nxm+j];
energy += sqrt(tmp.r*tmp.r+tmp.i*tmp.i);
}
}
}
if (verbose>1) fprintf(stderr,"eps_r=%e eps_a=%e\n", eps_r*energy, eps_a);
/* add small value at diagonal */
for (i=0; i<nxm; i++) {
tmp.r = eps_r*energy+eps_a;
matrix[i*nxm+i].r+=tmp.r;
}
jobz = "A";
lda = N = nxm;
lwork = N*N+4*N;
S = (float *)malloc(N*sizeof(float));
U = (complex *)malloc(N*N*sizeof(complex));
VT = (complex *)malloc(N*N*sizeof(complex));
work = (complex *)malloc(lwork*sizeof(complex));
rwork = (float *)malloc(5*(N*N+N)*sizeof(float));
iwork = (int *)malloc(8*N*sizeof(int));
/* Compute SVD */
cgesdd_(jobz, &N, &N, &matrix[0].r, &lda, S, &U[0].r, &lda, &VT[0].r,
&lda, &work[0].r, &lwork, rwork, iwork, &info);
assert (info == 0);
if (eigenvalues) {
for (i=0; i<N; i++) {
eigen[i] = S[i];
}
}
/* Compute inverse */
S0 = S[0];
neig = 0;
for (i=0; i<N; i++) {
/* fprintf(stderr,"S[%d] = %e S0 = %e\n ",i,S[i], S0);*/
Si = S[i];
if ((Si/S0) > numacc) { S[i]=1.0/S[i]; neig++; }
else S[i] = 0.0;
/* fprintf(stderr,"S^-1[%d] = %e\n",i,S[i]);*/
}
if(verbose) fprintf(stderr,"fraction of eigenvalues used = %.3f\n",(float)(neig/((float)N)));
for (j=0; j<N; j++) {
for (i=0; i<N; i++) {
U[j*N+i].r=S[j]*U[j*N+i].r;
U[j*N+i].i=-1.0*S[j]*U[j*N+i].i;
}
}
for (j=0; j<N; j++) {
for (i=0; i<N; i++) {
tmp.r = tmp.i = 0.0;
for (k=0; k<N; k++) {
a = U[k*N+j];
b.r = VT[i*N+k].r;
b.i = -1.0*VT[i*N+k].i;
tmp.r += (a.r*b.r-a.i*b.i);
tmp.i += (a.r*b.i+a.i*b.r);
}
matrix[j*nxm+i] = tmp;
}
}
free(U);
free(VT);
free(S);
free(work);
free(rwork);
free(iwork);
}
else if (rthm==4) { /* SVD general algorithm double precission most accurate */
double *rwork, *S, *U, *VT, ar, ai, br, bi, tmpr, tmpi;
double S0,Si,*Mat,*dwork;
int neig;
char *jobu, *jobvt;
energy=0.0;
if (eps_r != 0.0) {
for (i=0; i<nxm; i++) {
for (j=0; j<nxm; j++) {
tmp = matrix[i*nxm+j];
energy += sqrt(tmp.r*tmp.r+tmp.i*tmp.i);
}
}
}
if (verbose>1) fprintf(stderr,"eps_r=%e eps_a=%e\n", eps_r*energy, eps_a);
/* add small value at diagonal */
for (i=0; i<nxm; i++) {
tmp.r = eps_r*energy+eps_a;
matrix[i*nxm+i].r+=tmp.r;
}
Mat = (double *)malloc(2*N*N*sizeof(double));
/* convert to doubles */
for (i=0; i<nxm; i++) {
for (j=0; j<nxm; j++) {
Mat[i*2*nxm+j*2] = (double)matrix[i*nxm+j].r;
Mat[i*2*nxm+j*2+1] = (double)matrix[i*nxm+j].i;
}
}
jobu = "A";
jobvt = "A";
lda = N = nxm;
lwork = N*8;
S = (double *)malloc(N*sizeof(double));
U = (double *)malloc(2*N*N*sizeof(double));
VT = (double *)malloc(2*N*N*sizeof(double));
dwork = (double *)malloc(2*lwork*sizeof(double));
rwork = (double *)malloc(5*N*sizeof(double));
/* Compute SVD */
zgesvd_(jobu, jobvt, &N, &N, &Mat[0], &lda, S, &U[0], &lda, &VT[0],
&lda, &dwork[0], &lwork, rwork, &info);
assert (info == 0);
if (eigenvalues) {
for (i=0; i<N; i++) {
eigen[i] = (float)S[i];
}
}
/* Compute inverse */
S0 = S[0];
neig = 0;
for (i=0; i<N; i++) {
if (verbose=4) fprintf(stderr,"S[%d] = %e ",i,S[i]);
Si = S[i];
if ((Si/S0) > numacc) { S[i]=1.0/S[i]; neig++; }
else S[i] = 0.0;
/*S[i]=1.0/(S[i]+eps_r*S[0]);*/
/* fprintf(stderr,"S^-1[%d] = %e\n",i,S[i]);*/
}
if(verbose) fprintf(stderr,"fraction of eigenvalues used = %.3f\n",(float)(neig/((float)N)));
for (j=0; j<N; j++) {
for (i=0; i<N; i++) {
U[j*2*N+2*i]=S[j]*U[j*2*N+2*i];
U[j*2*N+2*i+1]=-1.0*S[j]*U[j*2*N+2*i+1];
}
}
for (j=0; j<N; j++) {
for (i=0; i<N; i++) {
tmpr = tmpi = 0.0;
for (k=0; k<N; k++) {
ar = U[k*2*N+2*j];
ai = U[k*2*N+2*j+1];
br = VT[i*2*N+2*k];
bi = -1.0*VT[i*2*N+2*k+1];
tmpr += (ar*br-ai*bi);
tmpi += (ar*bi+ai*br);
}
matrix[j*nxm+i].r = (float)tmpr;
matrix[j*nxm+i].i = (float)tmpi;
}
}
free(U);
free(VT);
free(S);
free(dwork);
free(rwork);
free(Mat);
}
else if (rthm==5) { /* double precission LU decomposition */
int ispec, n1, nb;
char *name , *opts;
double *Mat, *dwork;
ispec = 1;
name = "ZGETRI";
n1 = nxm;
nb = ilaenv_(&ispec, name, opts, &n1, &n1, &n1, &n1);
nb = MAX(1,nb);
lwork = nb*nxm;
ipvt = (int *)malloc(nxm*sizeof(int));
dwork = (double *)malloc(2*lwork*sizeof(double));
Mat = (double *)malloc(2*N*N*sizeof(double));
energy=0.0;
if (eps_r != 0.0) {
for (i=0; i<nxm; i++) {
for (j=0; j<nxm; j++) {
tmp = matrix[i*nxm+j];
energy += sqrt(tmp.r*tmp.r+tmp.i*tmp.i);
}
}
}
if (verbose>1) fprintf(stderr,"eps_r=%e eps_a=%e\n", eps_r*energy, eps_a);
/* convert to doubles */
for (i=0; i<nxm; i++) {
for (j=0; j<nxm; j++) {
Mat[i*2*nxm+j*2] = (double)matrix[i*nxm+j].r;
Mat[i*2*nxm+j*2+1] = (double)matrix[i*nxm+j].i;
}
}
/* add small value at diagonal */
for (i=0; i<nxm; i++) {
Mat[i*2*nxm+i*2] +=eps_r*energy+eps_a;
// Mat[i*2*nxm+i*2+1]+=eps_r*energy+eps_a;
}
/* LU based matrix inversion */
zgetrf_(&nxm, &nxm, &Mat[0], &nxm, ipvt, &info);
if (info != 0) fprintf(stderr,"error in zgetrf %d at frequency %d\n", info, iw);
assert (info == 0);
zgetri_(&nxm, &Mat[0], &nxm, ipvt, &dwork[0], &lwork, &info);
if (info != 0) fprintf(stderr,"error in zgetri %d at frequency %d\n", info, iw);
assert (info == 0);
/* convert back to floats */
for (i=0; i<nxm; i++) {
for (j=0; j<nxm; j++) {
matrix[i*nxm+j].r = (float)Mat[i*2*nxm+j*2];
matrix[i*nxm+j].i = (float)Mat[i*2*nxm+j*2+1];
}
}
free(ipvt);
free(dwork);
free(Mat);
}
else if (rthm==6) { /* eigenvalue decomposition */
int *iwork;
int neig;
double *work, *vr, *vl;
double *rwork, *S, *U, *VT, ar, ai, br, bi, tmpr, tmpi;
double S0,Si,nxi,*Mat;
char *jobvl, *jobvr;
jobvl = "V";
jobvr = "V";
lwork = N*N+2*N;
work = (double *)malloc(2*lwork*sizeof(double));
rwork = (double *)malloc(N*2*sizeof(double));
vr = (double *)malloc(2*N*N*sizeof(double));
vl = (double *)malloc(2*N*N*sizeof(double));
S = (double *)malloc(2*N*sizeof(double));
U = (double *)malloc(2*N*N*sizeof(double));
Mat = (double *)malloc(2*N*N*sizeof(double));
/* convert to doubles */
for (i=0; i<nxm; i++) {
for (j=0; j<nxm; j++) {
Mat[i*2*nxm+j*2] = (double)matrix[i*nxm+j].r;
Mat[i*2*nxm+j*2+1] = (double)matrix[i*nxm+j].i;
}
}
zgeev_(jobvl, jobvr, &N, Mat, &N, S, vl, &N, vr, &N,
work, &lwork, rwork, &info);
assert (info == 0);
nxi = 1.0/N;
for (i=0; i<N; i++) {
S[2*i] = (float)S[2*i]*nxi;
S[2*i+1] = (float)S[2*i+1]*nxi;
}
for (i=0; i<N; i++) {
for (j=0; j<N; j++) {
U[i*2*N+2*j] = (float)vr[(j)*2*N+2*i];
U[i*2*N+2*j+1] = (float)vr[(i)*2*N+2*j+1];
}
}
/* Compute inverse */
S0 = S[0];
neig = 0;
for (i=0; i<N; i++) {
/* fprintf(stderr,"S[%d] = %e ",i,S[i]);*/
Si = S[i];
if ((Si/S0) > numacc) { S[i]=1.0/S[i]; neig++; }
else S[i] = 0.0;
/* fprintf(stderr,"S^-1[%d] = %e\n",i,S[i]);*/
}
if(verbose) fprintf(stderr,"fraction of eigenvalues used = %.3f\n",(float)(neig/((float)N)));
for (j=0; j<N; j++) {
for (i=0; i<N; i++) {
U[j*2*N+2*i]=S[j]*U[j*2*N+2*i];
U[j*2*N+2*i+1]=-1.0*S[j]*U[j*2*N+2*i+1];
}
}
for (j=0; j<N; j++) {
for (i=0; i<N; i++) {
tmpr = tmpi = 0.0;
for (k=0; k<N; k++) {
ar = U[k*2*N+2*j];
ai = U[k*2*N+2*j+1];
br = U[i*2*N+2*k];
bi = U[i*2*N+2*k+1];
tmpr += (ar*br-ai*bi);
tmpi += (ar*bi+ai*br);
}
matrix[j*nxm+i].r = (float)tmpr;
matrix[j*nxm+i].i = (float)tmpi;
}
}
free(work);
free(rwork);
free(vr);
free(Mat);
free(S);
free(U);
}
return;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment