diff options
Diffstat (limited to 'src/lapack.cpp')
| -rw-r--r-- | src/lapack.cpp | 1030 |
1 files changed, 514 insertions, 516 deletions
diff --git a/src/lapack.cpp b/src/lapack.cpp index 05b85f4..8f6e8ff 100644 --- a/src/lapack.cpp +++ b/src/lapack.cpp @@ -16,614 +16,612 @@ along with this program. If not, see <http://www.gnu.org/licenses/>. */ -#include <iostream> +#include "gsl/gsl_linalg.h" +#include "gsl/gsl_matrix.h" +#include "gsl/gsl_vector.h" #include <cmath> +#include <iostream> #include <vector> -#include "gsl/gsl_vector.h" -#include "gsl/gsl_matrix.h" -#include "gsl/gsl_linalg.h" using namespace std; extern "C" void sgemm_(char *TRANSA, char *TRANSB, int *M, int *N, int *K, - float *ALPHA, float *A, int *LDA, float *B, int *LDB, - float *BETA, float *C, int *LDC); + float *ALPHA, float *A, int *LDA, float *B, int *LDB, + float *BETA, float *C, int *LDC); extern "C" void spotrf_(char *UPLO, int *N, float *A, int *LDA, int *INFO); extern "C" void spotrs_(char *UPLO, int *N, int *NRHS, float *A, int *LDA, - float *B, int *LDB, int *INFO); -extern "C" void ssyev_(char* JOBZ, char* UPLO, int *N, float *A, int *LDA, - float *W, float *WORK, int *LWORK, int *INFO); -extern "C" void ssyevr_(char* JOBZ, char *RANGE, char* UPLO, int *N, - float *A, int *LDA, float *VL, float *VU, int *IL, - int *IU, float *ABSTOL, int *M, float *W, float *Z, - int *LDZ, int *ISUPPZ, float *WORK, int *LWORK, - int *IWORK, int *LIWORK, int *INFO); + float *B, int *LDB, int *INFO); +extern "C" void ssyev_(char *JOBZ, char *UPLO, int *N, float *A, int *LDA, + float *W, float *WORK, int *LWORK, int *INFO); +extern "C" void ssyevr_(char *JOBZ, char *RANGE, char *UPLO, int *N, float *A, + int *LDA, float *VL, float *VU, int *IL, int *IU, + float *ABSTOL, int *M, float *W, float *Z, int *LDZ, + int *ISUPPZ, float *WORK, int *LWORK, int *IWORK, + int *LIWORK, int *INFO); extern "C" double sdot_(int *N, float *DX, int *INCX, float *DY, int *INCY); extern "C" void dgemm_(char *TRANSA, char *TRANSB, int *M, int *N, int *K, - double *ALPHA, double *A, int *LDA, double *B, - int *LDB, double *BETA, double *C, int *LDC); + double *ALPHA, double *A, int *LDA, double *B, int *LDB, + double *BETA, double *C, int *LDC); extern "C" void dpotrf_(char *UPLO, int *N, double *A, int *LDA, int *INFO); extern "C" void dpotrs_(char *UPLO, int *N, int *NRHS, double *A, int *LDA, - double *B, int *LDB, int *INFO); -extern "C" void dsyev_(char* JOBZ, char* UPLO, int *N, double *A, int *LDA, - double *W, double *WORK, int *LWORK, int *INFO); -extern "C" void dsyevr_(char* JOBZ, char *RANGE, char* UPLO, int *N, - double *A, int *LDA, double *VL, double *VU, - int *IL, int *IU, double *ABSTOL, int *M, - double *W, double *Z, int *LDZ, int *ISUPPZ, - double *WORK, int *LWORK, int *IWORK, - int *LIWORK, int *INFO); + double *B, int *LDB, int *INFO); +extern "C" void dsyev_(char *JOBZ, char *UPLO, int *N, double *A, int *LDA, + double *W, double *WORK, int *LWORK, int *INFO); +extern "C" void dsyevr_(char *JOBZ, char *RANGE, char *UPLO, int *N, double *A, + int *LDA, double *VL, double *VU, int *IL, int *IU, + double *ABSTOL, int *M, double *W, double *Z, int *LDZ, + int *ISUPPZ, double *WORK, int *LWORK, int *IWORK, + int *LIWORK, int *INFO); extern "C" double ddot_(int *N, double *DX, int *INCX, double *DY, int *INCY); // Cholesky decomposition, A is destroyed. -void lapack_float_cholesky_decomp (gsl_matrix_float *A) { - int N=A->size1, LDA=A->size1, INFO; - char UPLO='L'; - - if (N!=(int)A->size2) { - cout << "Matrix needs to be symmetric and same dimension in " << - "lapack_cholesky_decomp." << endl; - return; - } - - spotrf_(&UPLO, &N, A->data, &LDA, &INFO); - if (INFO!=0) { - cout << "Cholesky decomposition unsuccessful in " << - "lapack_cholesky_decomp." << endl; - return; - } - - return; +void lapack_float_cholesky_decomp(gsl_matrix_float *A) { + int N = A->size1, LDA = A->size1, INFO; + char UPLO = 'L'; + + if (N != (int)A->size2) { + cout << "Matrix needs to be symmetric and same dimension in " + << "lapack_cholesky_decomp." << endl; + return; + } + + spotrf_(&UPLO, &N, A->data, &LDA, &INFO); + if (INFO != 0) { + cout << "Cholesky decomposition unsuccessful in " + << "lapack_cholesky_decomp." << endl; + return; + } + + return; } // Cholesky decomposition, A is destroyed. -void lapack_cholesky_decomp (gsl_matrix *A) { - int N=A->size1, LDA=A->size1, INFO; - char UPLO='L'; - - if (N!=(int)A->size2) { - cout << "Matrix needs to be symmetric and same dimension in " << - "lapack_cholesky_decomp." << endl; - return; - } - - dpotrf_(&UPLO, &N, A->data, &LDA, &INFO); - if (INFO!=0) { - cout << "Cholesky decomposition unsuccessful in " << - "lapack_cholesky_decomp."<<endl; - return; - } - - return; +void lapack_cholesky_decomp(gsl_matrix *A) { + int N = A->size1, LDA = A->size1, INFO; + char UPLO = 'L'; + + if (N != (int)A->size2) { + cout << "Matrix needs to be symmetric and same dimension in " + << "lapack_cholesky_decomp." << endl; + return; + } + + dpotrf_(&UPLO, &N, A->data, &LDA, &INFO); + if (INFO != 0) { + cout << "Cholesky decomposition unsuccessful in " + << "lapack_cholesky_decomp." << endl; + return; + } + + return; } // Cholesky solve, A is decomposed. -void lapack_float_cholesky_solve (gsl_matrix_float *A, - const gsl_vector_float *b, - gsl_vector_float *x) { - int N=A->size1, NRHS=1, LDA=A->size1, LDB=b->size, INFO; - char UPLO='L'; - - - if (N!=(int)A->size2 || N!=LDB) { - cout << "Matrix needs to be symmetric and same dimension in " << - "lapack_cholesky_solve." << endl; - return; - } - - gsl_vector_float_memcpy (x, b); - spotrs_(&UPLO, &N, &NRHS, A->data, &LDA, x->data, &LDB, &INFO); - if (INFO!=0) { - cout << "Cholesky solve unsuccessful in lapack_cholesky_solve." << - endl; - return; - } - - return; +void lapack_float_cholesky_solve(gsl_matrix_float *A, const gsl_vector_float *b, + gsl_vector_float *x) { + int N = A->size1, NRHS = 1, LDA = A->size1, LDB = b->size, INFO; + char UPLO = 'L'; + + if (N != (int)A->size2 || N != LDB) { + cout << "Matrix needs to be symmetric and same dimension in " + << "lapack_cholesky_solve." << endl; + return; + } + + gsl_vector_float_memcpy(x, b); + spotrs_(&UPLO, &N, &NRHS, A->data, &LDA, x->data, &LDB, &INFO); + if (INFO != 0) { + cout << "Cholesky solve unsuccessful in lapack_cholesky_solve." << endl; + return; + } + + return; } // Cholesky solve, A is decomposed. -void lapack_cholesky_solve (gsl_matrix *A, const gsl_vector *b, - gsl_vector *x) { - int N=A->size1, NRHS=1, LDA=A->size1, LDB=b->size, INFO; - char UPLO='L'; - - if (N!=(int)A->size2 || N!=LDB) { - cout << "Matrix needs to be symmetric and same dimension in " << - "lapack_cholesky_solve." << endl; - return; - } - - gsl_vector_memcpy (x, b); - dpotrs_(&UPLO, &N, &NRHS, A->data, &LDA, x->data, &LDB, &INFO); - if (INFO!=0) { - cout << "Cholesky solve unsuccessful in lapack_cholesky_solve." << - endl; - return; - } - - return; -} +void lapack_cholesky_solve(gsl_matrix *A, const gsl_vector *b, gsl_vector *x) { + int N = A->size1, NRHS = 1, LDA = A->size1, LDB = b->size, INFO; + char UPLO = 'L'; + + if (N != (int)A->size2 || N != LDB) { + cout << "Matrix needs to be symmetric and same dimension in " + << "lapack_cholesky_solve." << endl; + return; + } + + gsl_vector_memcpy(x, b); + dpotrs_(&UPLO, &N, &NRHS, A->data, &LDA, x->data, &LDB, &INFO); + if (INFO != 0) { + cout << "Cholesky solve unsuccessful in lapack_cholesky_solve." << endl; + return; + } -void lapack_sgemm (char *TransA, char *TransB, float alpha, - const gsl_matrix_float *A, const gsl_matrix_float *B, - float beta, gsl_matrix_float *C) { - int M, N, K1, K2, LDA=A->size1, LDB=B->size1, LDC=C->size2; - - if (*TransA=='N' || *TransA=='n') {M=A->size1; K1=A->size2;} - else if (*TransA=='T' || *TransA=='t') {M=A->size2; K1=A->size1;} - else {cout<<"need 'N' or 'T' in lapack_sgemm"<<endl; return;} - - if (*TransB=='N' || *TransB=='n') {N=B->size2; K2=B->size1;} - else if (*TransB=='T' || *TransB=='t') {N=B->size1; K2=B->size2;} - else {cout<<"need 'N' or 'T' in lapack_sgemm"<<endl; return;} - - if (K1!=K2) { - cout<<"A and B not compatible in lapack_sgemm"<<endl; - return; - } - if (C->size1!=(size_t)M || C->size2!=(size_t)N) { - cout<<"C not compatible in lapack_sgemm"<<endl; - return; - } - - gsl_matrix_float *A_t=gsl_matrix_float_alloc (A->size2, A->size1); - gsl_matrix_float_transpose_memcpy (A_t, A); - gsl_matrix_float *B_t=gsl_matrix_float_alloc (B->size2, B->size1); - gsl_matrix_float_transpose_memcpy (B_t, B); - gsl_matrix_float *C_t=gsl_matrix_float_alloc (C->size2, C->size1); - gsl_matrix_float_transpose_memcpy (C_t, C); - - sgemm_(TransA, TransB, &M, &N, &K1, &alpha, A_t->data, &LDA, - B_t->data, &LDB, &beta, C_t->data, &LDC); - gsl_matrix_float_transpose_memcpy (C, C_t); - - gsl_matrix_float_free (A_t); - gsl_matrix_float_free (B_t); - gsl_matrix_float_free (C_t); - return; + return; } +void lapack_sgemm(char *TransA, char *TransB, float alpha, + const gsl_matrix_float *A, const gsl_matrix_float *B, + float beta, gsl_matrix_float *C) { + int M, N, K1, K2, LDA = A->size1, LDB = B->size1, LDC = C->size2; + + if (*TransA == 'N' || *TransA == 'n') { + M = A->size1; + K1 = A->size2; + } else if (*TransA == 'T' || *TransA == 't') { + M = A->size2; + K1 = A->size1; + } else { + cout << "need 'N' or 'T' in lapack_sgemm" << endl; + return; + } + if (*TransB == 'N' || *TransB == 'n') { + N = B->size2; + K2 = B->size1; + } else if (*TransB == 'T' || *TransB == 't') { + N = B->size1; + K2 = B->size2; + } else { + cout << "need 'N' or 'T' in lapack_sgemm" << endl; + return; + } -void lapack_dgemm (char *TransA, char *TransB, double alpha, - const gsl_matrix *A, const gsl_matrix *B, - double beta, gsl_matrix *C) { - int M, N, K1, K2, LDA=A->size1, LDB=B->size1, LDC=C->size2; + if (K1 != K2) { + cout << "A and B not compatible in lapack_sgemm" << endl; + return; + } + if (C->size1 != (size_t)M || C->size2 != (size_t)N) { + cout << "C not compatible in lapack_sgemm" << endl; + return; + } - if (*TransA=='N' || *TransA=='n') {M=A->size1; K1=A->size2;} - else if (*TransA=='T' || *TransA=='t') {M=A->size2; K1=A->size1;} - else {cout<<"need 'N' or 'T' in lapack_dgemm"<<endl; return;} + gsl_matrix_float *A_t = gsl_matrix_float_alloc(A->size2, A->size1); + gsl_matrix_float_transpose_memcpy(A_t, A); + gsl_matrix_float *B_t = gsl_matrix_float_alloc(B->size2, B->size1); + gsl_matrix_float_transpose_memcpy(B_t, B); + gsl_matrix_float *C_t = gsl_matrix_float_alloc(C->size2, C->size1); + gsl_matrix_float_transpose_memcpy(C_t, C); + + sgemm_(TransA, TransB, &M, &N, &K1, &alpha, A_t->data, &LDA, B_t->data, &LDB, + &beta, C_t->data, &LDC); + gsl_matrix_float_transpose_memcpy(C, C_t); + + gsl_matrix_float_free(A_t); + gsl_matrix_float_free(B_t); + gsl_matrix_float_free(C_t); + return; +} - if (*TransB=='N' || *TransB=='n') {N=B->size2; K2=B->size1;} - else if (*TransB=='T' || *TransB=='t') {N=B->size1; K2=B->size2;} - else {cout<<"need 'N' or 'T' in lapack_dgemm"<<endl; return;} +void lapack_dgemm(char *TransA, char *TransB, double alpha, const gsl_matrix *A, + const gsl_matrix *B, double beta, gsl_matrix *C) { + int M, N, K1, K2, LDA = A->size1, LDB = B->size1, LDC = C->size2; + + if (*TransA == 'N' || *TransA == 'n') { + M = A->size1; + K1 = A->size2; + } else if (*TransA == 'T' || *TransA == 't') { + M = A->size2; + K1 = A->size1; + } else { + cout << "need 'N' or 'T' in lapack_dgemm" << endl; + return; + } + + if (*TransB == 'N' || *TransB == 'n') { + N = B->size2; + K2 = B->size1; + } else if (*TransB == 'T' || *TransB == 't') { + N = B->size1; + K2 = B->size2; + } else { + cout << "need 'N' or 'T' in lapack_dgemm" << endl; + return; + } - if (K1!=K2) { - cout << "A and B not compatible in lapack_dgemm"<<endl; - return; - } - if (C->size1!=(size_t)M || C->size2!=(size_t)N) { - cout<<"C not compatible in lapack_dgemm"<<endl; - return; - } + if (K1 != K2) { + cout << "A and B not compatible in lapack_dgemm" << endl; + return; + } + if (C->size1 != (size_t)M || C->size2 != (size_t)N) { + cout << "C not compatible in lapack_dgemm" << endl; + return; + } - gsl_matrix *A_t=gsl_matrix_alloc (A->size2, A->size1); - gsl_matrix_transpose_memcpy (A_t, A); - gsl_matrix *B_t=gsl_matrix_alloc (B->size2, B->size1); - gsl_matrix_transpose_memcpy (B_t, B); - gsl_matrix *C_t=gsl_matrix_alloc (C->size2, C->size1); - gsl_matrix_transpose_memcpy (C_t, C); + gsl_matrix *A_t = gsl_matrix_alloc(A->size2, A->size1); + gsl_matrix_transpose_memcpy(A_t, A); + gsl_matrix *B_t = gsl_matrix_alloc(B->size2, B->size1); + gsl_matrix_transpose_memcpy(B_t, B); + gsl_matrix *C_t = gsl_matrix_alloc(C->size2, C->size1); + gsl_matrix_transpose_memcpy(C_t, C); - dgemm_(TransA, TransB, &M, &N, &K1, &alpha, A_t->data, &LDA, - B_t->data, &LDB, &beta, C_t->data, &LDC); + dgemm_(TransA, TransB, &M, &N, &K1, &alpha, A_t->data, &LDA, B_t->data, &LDB, + &beta, C_t->data, &LDC); - gsl_matrix_transpose_memcpy (C, C_t); + gsl_matrix_transpose_memcpy(C, C_t); - gsl_matrix_free (A_t); - gsl_matrix_free (B_t); - gsl_matrix_free (C_t); - return; + gsl_matrix_free(A_t); + gsl_matrix_free(B_t); + gsl_matrix_free(C_t); + return; } // Eigen value decomposition, matrix A is destroyed, float seems to // have problem with large matrices (in mac). -void lapack_float_eigen_symmv (gsl_matrix_float *A, gsl_vector_float *eval, - gsl_matrix_float *evec, - const size_t flag_largematrix) { - if (flag_largematrix==1) { - int N=A->size1, LDA=A->size1, INFO, LWORK=-1; - char JOBZ='V', UPLO='L'; - - if (N!=(int)A->size2 || N!=(int)eval->size) { - cout << "Matrix needs to be symmetric and same " << - "dimension in lapack_eigen_symmv."<<endl; - return; - } - - LWORK=3*N; - float *WORK=new float [LWORK]; - ssyev_(&JOBZ, &UPLO, &N, A->data, &LDA, eval->data, WORK, - &LWORK, &INFO); - if (INFO!=0) { - cout << "Eigen decomposition unsuccessful in " << - "lapack_eigen_symmv."<<endl; - return; - } - - gsl_matrix_float_view A_sub = - gsl_matrix_float_submatrix(A, 0, 0, N, N); - gsl_matrix_float_memcpy (evec, &A_sub.matrix); - gsl_matrix_float_transpose (evec); - - delete [] WORK; - } else { - int N=A->size1, LDA=A->size1, LDZ=A->size1, INFO, - LWORK=-1, LIWORK=-1; - char JOBZ='V', UPLO='L', RANGE='A'; - float ABSTOL=1.0E-7; - - // VL, VU, IL, IU are not referenced; M equals N if RANGE='A'. - float VL=0.0, VU=0.0; - int IL=0, IU=0, M; - - if (N!=(int)A->size2 || N!=(int)eval->size) { - cout << "Matrix needs to be symmetric and same " << - "dimension in lapack_float_eigen_symmv." << endl; - return; - } - - int *ISUPPZ=new int [2*N]; - - float WORK_temp[1]; - int IWORK_temp[1]; - ssyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, - &VU, &IL, &IU, &ABSTOL, &M, eval->data, - evec->data, &LDZ, ISUPPZ, WORK_temp, &LWORK, - IWORK_temp, &LIWORK, &INFO); - if (INFO!=0) { - cout << "Work space estimate unsuccessful in " << - "lapack_float_eigen_symmv." << endl; - return; - } - LWORK=(int)WORK_temp[0]; LIWORK=(int)IWORK_temp[0]; - - float *WORK=new float [LWORK]; - int *IWORK=new int [LIWORK]; - - ssyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, - &VU, &IL, &IU, &ABSTOL, &M, eval->data, evec->data, - &LDZ, ISUPPZ, WORK, &LWORK, IWORK, &LIWORK, &INFO); - if (INFO!=0) { - cout << "Eigen decomposition unsuccessful in " << - "lapack_float_eigen_symmv." << endl; - return; - } - - gsl_matrix_float_transpose (evec); - - delete [] ISUPPZ; - delete [] WORK; - delete [] IWORK; - } - - - return; -} - +void lapack_float_eigen_symmv(gsl_matrix_float *A, gsl_vector_float *eval, + gsl_matrix_float *evec, + const size_t flag_largematrix) { + if (flag_largematrix == 1) { + int N = A->size1, LDA = A->size1, INFO, LWORK = -1; + char JOBZ = 'V', UPLO = 'L'; + + if (N != (int)A->size2 || N != (int)eval->size) { + cout << "Matrix needs to be symmetric and same " + << "dimension in lapack_eigen_symmv." << endl; + return; + } + + LWORK = 3 * N; + float *WORK = new float[LWORK]; + ssyev_(&JOBZ, &UPLO, &N, A->data, &LDA, eval->data, WORK, &LWORK, &INFO); + if (INFO != 0) { + cout << "Eigen decomposition unsuccessful in " + << "lapack_eigen_symmv." << endl; + return; + } + + gsl_matrix_float_view A_sub = gsl_matrix_float_submatrix(A, 0, 0, N, N); + gsl_matrix_float_memcpy(evec, &A_sub.matrix); + gsl_matrix_float_transpose(evec); + + delete[] WORK; + } else { + int N = A->size1, LDA = A->size1, LDZ = A->size1, INFO, LWORK = -1, + LIWORK = -1; + char JOBZ = 'V', UPLO = 'L', RANGE = 'A'; + float ABSTOL = 1.0E-7; + + // VL, VU, IL, IU are not referenced; M equals N if RANGE='A'. + float VL = 0.0, VU = 0.0; + int IL = 0, IU = 0, M; + + if (N != (int)A->size2 || N != (int)eval->size) { + cout << "Matrix needs to be symmetric and same " + << "dimension in lapack_float_eigen_symmv." << endl; + return; + } + + int *ISUPPZ = new int[2 * N]; + + float WORK_temp[1]; + int IWORK_temp[1]; + ssyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, &VU, &IL, &IU, + &ABSTOL, &M, eval->data, evec->data, &LDZ, ISUPPZ, WORK_temp, + &LWORK, IWORK_temp, &LIWORK, &INFO); + if (INFO != 0) { + cout << "Work space estimate unsuccessful in " + << "lapack_float_eigen_symmv." << endl; + return; + } + LWORK = (int)WORK_temp[0]; + LIWORK = (int)IWORK_temp[0]; + + float *WORK = new float[LWORK]; + int *IWORK = new int[LIWORK]; + + ssyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, &VU, &IL, &IU, + &ABSTOL, &M, eval->data, evec->data, &LDZ, ISUPPZ, WORK, &LWORK, + IWORK, &LIWORK, &INFO); + if (INFO != 0) { + cout << "Eigen decomposition unsuccessful in " + << "lapack_float_eigen_symmv." << endl; + return; + } + + gsl_matrix_float_transpose(evec); + + delete[] ISUPPZ; + delete[] WORK; + delete[] IWORK; + } + return; +} // Eigenvalue decomposition, matrix A is destroyed. -void lapack_eigen_symmv (gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec, - const size_t flag_largematrix) { - if (flag_largematrix==1) { - int N=A->size1, LDA=A->size1, INFO, LWORK=-1; - char JOBZ='V', UPLO='L'; - - if (N!=(int)A->size2 || N!=(int)eval->size) { - cout << "Matrix needs to be symmetric and same " << - "dimension in lapack_eigen_symmv." << endl; - return; - } - - LWORK=3*N; - double *WORK=new double [LWORK]; - dsyev_(&JOBZ, &UPLO, &N, A->data, &LDA, eval->data, WORK, - &LWORK, &INFO); - if (INFO!=0) { - cout<<"Eigen decomposition unsuccessful in " << - "lapack_eigen_symmv." << endl; - return; - } - - gsl_matrix_view A_sub=gsl_matrix_submatrix(A, 0, 0, N, N); - gsl_matrix_memcpy (evec, &A_sub.matrix); - gsl_matrix_transpose (evec); - - delete [] WORK; - } else { - int N=A->size1, LDA=A->size1, LDZ=A->size1, INFO; - int LWORK=-1, LIWORK=-1; - char JOBZ='V', UPLO='L', RANGE='A'; - double ABSTOL=1.0E-7; - - // VL, VU, IL, IU are not referenced; M equals N if RANGE='A'. - double VL=0.0, VU=0.0; - int IL=0, IU=0, M; - - if (N!=(int)A->size2 || N!=(int)eval->size) { - cout << "Matrix needs to be symmetric and same " << - "dimension in lapack_eigen_symmv." << endl; - return; - } - - int *ISUPPZ=new int [2*N]; - - double WORK_temp[1]; - int IWORK_temp[1]; - - dsyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, &VU, - &IL, &IU, &ABSTOL, &M, eval->data, evec->data, - &LDZ, ISUPPZ, WORK_temp, &LWORK, IWORK_temp, - &LIWORK, &INFO); - if (INFO!=0) { - cout << "Work space estimate unsuccessful in " << - "lapack_eigen_symmv." << endl; - return; - } - LWORK=(int)WORK_temp[0]; LIWORK=(int)IWORK_temp[0]; - - double *WORK=new double [LWORK]; - int *IWORK=new int [LIWORK]; - - dsyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, &VU, - &IL, &IU, &ABSTOL, &M, eval->data, evec->data, - &LDZ, ISUPPZ, WORK, &LWORK, IWORK, &LIWORK, &INFO); - if (INFO!=0) { - cout << "Eigen decomposition unsuccessful in " << - "lapack_eigen_symmv." << endl; - return; - } - - gsl_matrix_transpose (evec); - - delete [] ISUPPZ; - delete [] WORK; - delete [] IWORK; - } - - return; +void lapack_eigen_symmv(gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec, + const size_t flag_largematrix) { + if (flag_largematrix == 1) { + int N = A->size1, LDA = A->size1, INFO, LWORK = -1; + char JOBZ = 'V', UPLO = 'L'; + + if (N != (int)A->size2 || N != (int)eval->size) { + cout << "Matrix needs to be symmetric and same " + << "dimension in lapack_eigen_symmv." << endl; + return; + } + + LWORK = 3 * N; + double *WORK = new double[LWORK]; + dsyev_(&JOBZ, &UPLO, &N, A->data, &LDA, eval->data, WORK, &LWORK, &INFO); + if (INFO != 0) { + cout << "Eigen decomposition unsuccessful in " + << "lapack_eigen_symmv." << endl; + return; + } + + gsl_matrix_view A_sub = gsl_matrix_submatrix(A, 0, 0, N, N); + gsl_matrix_memcpy(evec, &A_sub.matrix); + gsl_matrix_transpose(evec); + + delete[] WORK; + } else { + int N = A->size1, LDA = A->size1, LDZ = A->size1, INFO; + int LWORK = -1, LIWORK = -1; + char JOBZ = 'V', UPLO = 'L', RANGE = 'A'; + double ABSTOL = 1.0E-7; + + // VL, VU, IL, IU are not referenced; M equals N if RANGE='A'. + double VL = 0.0, VU = 0.0; + int IL = 0, IU = 0, M; + + if (N != (int)A->size2 || N != (int)eval->size) { + cout << "Matrix needs to be symmetric and same " + << "dimension in lapack_eigen_symmv." << endl; + return; + } + + int *ISUPPZ = new int[2 * N]; + + double WORK_temp[1]; + int IWORK_temp[1]; + + dsyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, &VU, &IL, &IU, + &ABSTOL, &M, eval->data, evec->data, &LDZ, ISUPPZ, WORK_temp, + &LWORK, IWORK_temp, &LIWORK, &INFO); + if (INFO != 0) { + cout << "Work space estimate unsuccessful in " + << "lapack_eigen_symmv." << endl; + return; + } + LWORK = (int)WORK_temp[0]; + LIWORK = (int)IWORK_temp[0]; + + double *WORK = new double[LWORK]; + int *IWORK = new int[LIWORK]; + + dsyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, &VU, &IL, &IU, + &ABSTOL, &M, eval->data, evec->data, &LDZ, ISUPPZ, WORK, &LWORK, + IWORK, &LIWORK, &INFO); + if (INFO != 0) { + cout << "Eigen decomposition unsuccessful in " + << "lapack_eigen_symmv." << endl; + return; + } + + gsl_matrix_transpose(evec); + + delete[] ISUPPZ; + delete[] WORK; + delete[] IWORK; + } + + return; } // DO NOT set eigenvalues to be positive. -double EigenDecomp (gsl_matrix *G, gsl_matrix *U, gsl_vector *eval, - const size_t flag_largematrix) { - lapack_eigen_symmv (G, eval, U, flag_largematrix); - - // Calculate track_G=mean(diag(G)). - double d=0.0; - for (size_t i=0; i<eval->size; ++i) { - d+=gsl_vector_get(eval, i); - } - d/=(double)eval->size; - - return d; -} +double EigenDecomp(gsl_matrix *G, gsl_matrix *U, gsl_vector *eval, + const size_t flag_largematrix) { + lapack_eigen_symmv(G, eval, U, flag_largematrix); + + // Calculate track_G=mean(diag(G)). + double d = 0.0; + for (size_t i = 0; i < eval->size; ++i) { + d += gsl_vector_get(eval, i); + } + d /= (double)eval->size; + return d; +} // DO NOT set eigen values to be positive. -double EigenDecomp (gsl_matrix_float *G, gsl_matrix_float *U, - gsl_vector_float *eval, const size_t flag_largematrix) { - lapack_float_eigen_symmv (G, eval, U, flag_largematrix); - - // Calculate track_G=mean(diag(G)). - double d = 0.0; - for (size_t i=0; i<eval->size; ++i) { - d+=gsl_vector_float_get(eval, i); - } - d/=(double)eval->size; - - return d; -} +double EigenDecomp(gsl_matrix_float *G, gsl_matrix_float *U, + gsl_vector_float *eval, const size_t flag_largematrix) { + lapack_float_eigen_symmv(G, eval, U, flag_largematrix); + + // Calculate track_G=mean(diag(G)). + double d = 0.0; + for (size_t i = 0; i < eval->size; ++i) { + d += gsl_vector_float_get(eval, i); + } + d /= (double)eval->size; + return d; +} double CholeskySolve(gsl_matrix *Omega, gsl_vector *Xty, gsl_vector *OiXty) { - double logdet_O=0.0; + double logdet_O = 0.0; - lapack_cholesky_decomp(Omega); - for (size_t i=0; i<Omega->size1; ++i) { - logdet_O+=log(gsl_matrix_get (Omega, i, i)); - } - logdet_O*=2.0; - lapack_cholesky_solve(Omega, Xty, OiXty); + lapack_cholesky_decomp(Omega); + for (size_t i = 0; i < Omega->size1; ++i) { + logdet_O += log(gsl_matrix_get(Omega, i, i)); + } + logdet_O *= 2.0; + lapack_cholesky_solve(Omega, Xty, OiXty); - return logdet_O; + return logdet_O; } - double CholeskySolve(gsl_matrix_float *Omega, gsl_vector_float *Xty, - gsl_vector_float *OiXty) { - double logdet_O=0.0; + gsl_vector_float *OiXty) { + double logdet_O = 0.0; - lapack_float_cholesky_decomp(Omega); - for (size_t i=0; i<Omega->size1; ++i) { - logdet_O+=log(gsl_matrix_float_get (Omega, i, i)); - } - logdet_O*=2.0; - lapack_float_cholesky_solve(Omega, Xty, OiXty); + lapack_float_cholesky_decomp(Omega); + for (size_t i = 0; i < Omega->size1; ++i) { + logdet_O += log(gsl_matrix_float_get(Omega, i, i)); + } + logdet_O *= 2.0; + lapack_float_cholesky_solve(Omega, Xty, OiXty); - return logdet_O; + return logdet_O; } - // LU decomposition. -void LUDecomp (gsl_matrix *LU, gsl_permutation *p, int *signum) { - gsl_linalg_LU_decomp (LU, p, signum); - return; +void LUDecomp(gsl_matrix *LU, gsl_permutation *p, int *signum) { + gsl_linalg_LU_decomp(LU, p, signum); + return; } -void LUDecomp (gsl_matrix_float *LU, gsl_permutation *p, int *signum) { - gsl_matrix *LU_double=gsl_matrix_alloc (LU->size1, LU->size2); - - // Copy float matrix to double. - for (size_t i=0; i<LU->size1; i++) { - for (size_t j=0; j<LU->size2; j++) { - gsl_matrix_set (LU_double, i, j, - gsl_matrix_float_get(LU, i, j)); - } - } - - // LU decomposition. - gsl_linalg_LU_decomp (LU_double, p, signum); - - // Copy float matrix to double. - for (size_t i=0; i<LU->size1; i++) { - for (size_t j=0; j<LU->size2; j++) { - gsl_matrix_float_set (LU, i, j, - gsl_matrix_get(LU_double, i, j)); - } - } - - // Free matrix. - gsl_matrix_free (LU_double); - return; -} +void LUDecomp(gsl_matrix_float *LU, gsl_permutation *p, int *signum) { + gsl_matrix *LU_double = gsl_matrix_alloc(LU->size1, LU->size2); + + // Copy float matrix to double. + for (size_t i = 0; i < LU->size1; i++) { + for (size_t j = 0; j < LU->size2; j++) { + gsl_matrix_set(LU_double, i, j, gsl_matrix_float_get(LU, i, j)); + } + } + // LU decomposition. + gsl_linalg_LU_decomp(LU_double, p, signum); + + // Copy float matrix to double. + for (size_t i = 0; i < LU->size1; i++) { + for (size_t j = 0; j < LU->size2; j++) { + gsl_matrix_float_set(LU, i, j, gsl_matrix_get(LU_double, i, j)); + } + } + + // Free matrix. + gsl_matrix_free(LU_double); + return; +} // LU invert. -void LUInvert (const gsl_matrix *LU, const gsl_permutation *p, - gsl_matrix *inverse) { - gsl_linalg_LU_invert (LU, p, inverse); - return; +void LUInvert(const gsl_matrix *LU, const gsl_permutation *p, + gsl_matrix *inverse) { + gsl_linalg_LU_invert(LU, p, inverse); + return; } -void LUInvert (const gsl_matrix_float *LU, const gsl_permutation *p, - gsl_matrix_float *inverse) { - gsl_matrix *LU_double=gsl_matrix_alloc (LU->size1, LU->size2); - gsl_matrix *inverse_double=gsl_matrix_alloc (inverse->size1, - inverse->size2); - - // Copy float matrix to double. - for (size_t i=0; i<LU->size1; i++) { - for (size_t j=0; j<LU->size2; j++) { - gsl_matrix_set (LU_double, i, j, - gsl_matrix_float_get(LU, i, j)); - } - } - - // LU decomposition. - gsl_linalg_LU_invert (LU_double, p, inverse_double); - - // Copy float matrix to double. - for (size_t i=0; i<inverse->size1; i++) { - for (size_t j=0; j<inverse->size2; j++) { - gsl_matrix_float_set (inverse, i, j, - gsl_matrix_get(inverse_double, - i, j)); - } - } - - // Free matrix. - gsl_matrix_free (LU_double); - gsl_matrix_free (inverse_double); - return; +void LUInvert(const gsl_matrix_float *LU, const gsl_permutation *p, + gsl_matrix_float *inverse) { + gsl_matrix *LU_double = gsl_matrix_alloc(LU->size1, LU->size2); + gsl_matrix *inverse_double = gsl_matrix_alloc(inverse->size1, inverse->size2); + + // Copy float matrix to double. + for (size_t i = 0; i < LU->size1; i++) { + for (size_t j = 0; j < LU->size2; j++) { + gsl_matrix_set(LU_double, i, j, gsl_matrix_float_get(LU, i, j)); + } + } + + // LU decomposition. + gsl_linalg_LU_invert(LU_double, p, inverse_double); + + // Copy float matrix to double. + for (size_t i = 0; i < inverse->size1; i++) { + for (size_t j = 0; j < inverse->size2; j++) { + gsl_matrix_float_set(inverse, i, j, gsl_matrix_get(inverse_double, i, j)); + } + } + + // Free matrix. + gsl_matrix_free(LU_double); + gsl_matrix_free(inverse_double); + return; } // LU lndet. -double LULndet (gsl_matrix *LU) { - double d; - d=gsl_linalg_LU_lndet (LU); - return d; +double LULndet(gsl_matrix *LU) { + double d; + d = gsl_linalg_LU_lndet(LU); + return d; } -double LULndet (gsl_matrix_float *LU) { - gsl_matrix *LU_double=gsl_matrix_alloc (LU->size1, LU->size2); - double d; +double LULndet(gsl_matrix_float *LU) { + gsl_matrix *LU_double = gsl_matrix_alloc(LU->size1, LU->size2); + double d; - // Copy float matrix to double. - for (size_t i=0; i<LU->size1; i++) { - for (size_t j=0; j<LU->size2; j++) { - gsl_matrix_set (LU_double, i, j, gsl_matrix_float_get(LU, i, j)); - } - } + // Copy float matrix to double. + for (size_t i = 0; i < LU->size1; i++) { + for (size_t j = 0; j < LU->size2; j++) { + gsl_matrix_set(LU_double, i, j, gsl_matrix_float_get(LU, i, j)); + } + } - // LU decomposition. - d=gsl_linalg_LU_lndet (LU_double); + // LU decomposition. + d = gsl_linalg_LU_lndet(LU_double); - // Free matrix - gsl_matrix_free (LU_double); - return d; + // Free matrix + gsl_matrix_free(LU_double); + return d; } - // LU solve. -void LUSolve (const gsl_matrix *LU, const gsl_permutation *p, - const gsl_vector *b, gsl_vector *x) { - gsl_linalg_LU_solve (LU, p, b, x); - return; +void LUSolve(const gsl_matrix *LU, const gsl_permutation *p, + const gsl_vector *b, gsl_vector *x) { + gsl_linalg_LU_solve(LU, p, b, x); + return; } -void LUSolve (const gsl_matrix_float *LU, const gsl_permutation *p, - const gsl_vector_float *b, gsl_vector_float *x) { - gsl_matrix *LU_double=gsl_matrix_alloc (LU->size1, LU->size2); - gsl_vector *b_double=gsl_vector_alloc (b->size); - gsl_vector *x_double=gsl_vector_alloc (x->size); - - // Copy float matrix to double. - for (size_t i=0; i<LU->size1; i++) { - for (size_t j=0; j<LU->size2; j++) { - gsl_matrix_set (LU_double, i, j, - gsl_matrix_float_get(LU, i, j)); - } - } - - for (size_t i=0; i<b->size; i++) { - gsl_vector_set (b_double, i, gsl_vector_float_get(b, i)); - } - - for (size_t i=0; i<x->size; i++) { - gsl_vector_set (x_double, i, gsl_vector_float_get(x, i)); - } - - // LU decomposition. - gsl_linalg_LU_solve (LU_double, p, b_double, x_double); - - // Copy float matrix to double. - for (size_t i=0; i<x->size; i++) { - gsl_vector_float_set (x, i, gsl_vector_get(x_double, i)); - } - - // Free matrix. - gsl_matrix_free (LU_double); - gsl_vector_free (b_double); - gsl_vector_free (x_double); - return; -} +void LUSolve(const gsl_matrix_float *LU, const gsl_permutation *p, + const gsl_vector_float *b, gsl_vector_float *x) { + gsl_matrix *LU_double = gsl_matrix_alloc(LU->size1, LU->size2); + gsl_vector *b_double = gsl_vector_alloc(b->size); + gsl_vector *x_double = gsl_vector_alloc(x->size); + + // Copy float matrix to double. + for (size_t i = 0; i < LU->size1; i++) { + for (size_t j = 0; j < LU->size2; j++) { + gsl_matrix_set(LU_double, i, j, gsl_matrix_float_get(LU, i, j)); + } + } + for (size_t i = 0; i < b->size; i++) { + gsl_vector_set(b_double, i, gsl_vector_float_get(b, i)); + } + + for (size_t i = 0; i < x->size; i++) { + gsl_vector_set(x_double, i, gsl_vector_float_get(x, i)); + } + + // LU decomposition. + gsl_linalg_LU_solve(LU_double, p, b_double, x_double); + + // Copy float matrix to double. + for (size_t i = 0; i < x->size; i++) { + gsl_vector_float_set(x, i, gsl_vector_get(x_double, i)); + } + + // Free matrix. + gsl_matrix_free(LU_double); + gsl_vector_free(b_double); + gsl_vector_free(x_double); + return; +} bool lapack_ddot(vector<double> &x, vector<double> &y, double &v) { - bool flag=false; - int incx=1; - int incy=1; - int n=(int)x.size(); - if (x.size()==y.size()) { - v=ddot_(&n, &x[0], &incx, &y[0], &incy); - flag=true; + bool flag = false; + int incx = 1; + int incy = 1; + int n = (int)x.size(); + if (x.size() == y.size()) { + v = ddot_(&n, &x[0], &incx, &y[0], &incy); + flag = true; } return flag; } - bool lapack_sdot(vector<float> &x, vector<float> &y, double &v) { - bool flag=false; - int incx=1; - int incy=1; - int n=(int)x.size(); - if (x.size()==y.size()) { - v=sdot_(&n, &x[0], &incx, &y[0], &incy); - flag=true; + bool flag = false; + int incx = 1; + int incy = 1; + int n = (int)x.size(); + if (x.size() == y.size()) { + v = sdot_(&n, &x[0], &incx, &y[0], &incy); + flag = true; } return flag; |