1 /* Interface with the LAPACK (Fortran77) numerical library.
2 */
3 #include "esl_config.h"
4 #ifdef HAVE_LIBLAPACK
5
6 #include <stdlib.h>
7
8 #include "easel.h"
9 #include "esl_dmatrix.h"
10 #include "interface_lapack.h"
11
12 /* A: nxn real matrix
13 * ret_Er: RETURN: vector of eigenvalues, real part, allocated 0..n-1
14 * ret_Ei: RETURN: vector of eigenvalues, imaginary part, allocated 0..n-1
15 * ret_VL: RETURN: left eigenvectors
16 * ret_VR: RETURN: right eigenvectors
17 */
18 int
esl_lapack_dgeev(ESL_DMATRIX * A,double ** ret_Er,double ** ret_Ei,ESL_DMATRIX ** ret_VL,ESL_DMATRIX ** ret_VR)19 esl_lapack_dgeev(ESL_DMATRIX *A, double **ret_Er, double **ret_Ei, ESL_DMATRIX **ret_VL, ESL_DMATRIX **ret_VR)
20 {
21 double *Er = NULL;
22 double *Ei = NULL;
23 ESL_DMATRIX *VL = NULL;
24 ESL_DMATRIX *VR = NULL;
25 double *work = NULL;
26 char jobvl, jobvr;
27 int lda;
28 int ldvl, ldvr;
29 int lwork;
30 int info;
31 int status;
32
33 if ((VL = esl_dmatrix_Create(A->n,A->n)) == NULL) { status = eslEMEM; goto ERROR; }
34 if ((VR = esl_dmatrix_Create(A->n,A->n)) == NULL) { status = eslEMEM; goto ERROR; }
35 ESL_ALLOC(Er, sizeof(double) * A->n);
36 ESL_ALLOC(Ei, sizeof(double) * A->n);
37 ESL_ALLOC(work, sizeof(double) * 4 * A->n);
38
39 jobvl = (ret_VL == NULL) ? 'N' : 'V'; /* do we want left eigenvectors? */
40 jobvr = (ret_VR == NULL) ? 'N' : 'V'; /* do we want right eigenvectors? */
41 lda = A->n;
42 ldvl = A->n;
43 ldvr = A->n;
44 lwork = 4*A->n;
45
46 /* Fortran convention is colxrow, not rowxcol; so transpose
47 * A before passing it to a Fortran routine.
48 */
49 esl_dmx_Transpose(A);
50
51 /* The actual Fortran77 interface call to LAPACK.
52 * All args must be passed by reference.
53 * Fortran 2D arrays are 1D: so pass the A[0] part of a DSMX.
54 */
55 dgeev_(&jobvl, &jobvr, &(A->n), A->mx[0], &lda, Er, Ei, VL->mx[0], &ldvl, VR->mx[0], &ldvr, work, &lwork, &info);
56
57 /* Now, VL, VR are transposed (col x row), so transpose them back to
58 * C convention.
59 */
60 esl_dmx_Transpose(VL);
61 esl_dmx_Transpose(VR);
62
63 if (ret_VL != NULL) *ret_VL = VL; else esl_dmatrix_Destroy(VL);
64 if (ret_VR != NULL) *ret_VR = VR; else esl_dmatrix_Destroy(VR);
65 if (ret_Er != NULL) *ret_Er = Er; else free(Er);
66 if (ret_Ei != NULL) *ret_Ei = Ei; else free(Ei);
67 free(work);
68 return eslOK;
69
70 ERROR:
71 if (ret_VL != NULL) *ret_VL = NULL;
72 if (ret_VR != NULL) *ret_VR = NULL;
73 if (ret_Er != NULL) *ret_Er = NULL;
74 if (ret_Ei != NULL) *ret_Ei = NULL;
75 if (VL != NULL) free(VL);
76 if (VR != NULL) free(VR);
77 if (Er != NULL) free(Er);
78 if (Ei != NULL) free(Ei);
79 if (work != NULL) free(work);
80 return status;
81 }
82
83
84 #endif /*HAVE_LIBLAPACK*/
85
86
87