1 /* _______________________________________________________________________ 2 3 Surfpack: A Software Library of Multidimensional Surface Fitting Methods 4 Copyright (c) 2006, Sandia National Laboratories. 5 This software is distributed under the GNU Lesser General Public License. 6 For more information, see the README file in the top Surfpack directory. 7 _______________________________________________________________________ */ 8 9 #ifndef __SURFPACK_LAPACK_WRAPPERS_H__ 10 #define __SURFPACK_LAPACK_WRAPPERS_H__ 11 12 #ifdef HAVE_CONFIG_H 13 #include "surfpack_config.h" 14 #else 15 #include "surf77_config.h" 16 #endif 17 18 /***************************************************************************/ 19 /***************************************************************************/ 20 /**** The BLAS and LAPACK wrappers should be the same whichever version ****/ 21 /**** of SurfMat we use, and since they are _ONLY_ wrappers they should ****/ 22 /**** the should be inline functions so they should be in a header file ****/ 23 /**** so it makes sense to put them here. ****/ 24 /***************************************************************************/ 25 /***************************************************************************/ 26 27 #ifdef HAVE_CONFIG_H 28 29 // BMA (2016): Changed to use Fortran 2003 ISO C bindings. 30 // The Fortran symbol will be lowercase with same name as if in C 31 //#define NKM_PIVOTCHOL_F77 F77_FUNC(nkm_pivotchol,NKM_PIVOTCHOL) 32 //#define NKM_BLOCKPIVOTCHOL_F77 F77_FUNC(nkm_blockpivotchol,NKM_BLOCKPIVOTCHOL) 33 #define NKM_PIVOTCHOL_F77 nkm_pivotchol 34 #define NKM_BLOCKPIVOTCHOL_F77 nkm_blockpivotchol 35 36 37 #define DGETRF_F77 F77_FUNC(dgetrf,DGETRF) 38 #define DGETRI_F77 F77_FUNC(dgetri,DGETRI) 39 #define DGEMV_F77 F77_FUNC(dgemv,DGEMV) 40 #define DGEMM_F77 F77_FUNC(dgemm,DGEMM) 41 #define DDOT_F77 F77_FUNC(ddot, DDOT) 42 #define DGELS_F77 F77_FUNC(dgels,DGELS) 43 #define DGESVD_F77 F77_FUNC(dgesvd,DGESVD) 44 45 #define DPOTRF_F77 F77_FUNC(dpotrf,DPOTRF) 46 #define DPOTRI_F77 F77_FUNC(dpotri,DPOTRI) 47 #define DPOTRS_F77 F77_FUNC(dpotrs,DPOTRS) 48 #define DPOCON_F77 F77_FUNC(dpocon,DPOCON) 49 #define DGETRS_F77 F77_FUNC(dgetrs,DGETRS) 50 #define DLANGE_F77 F77_FUNC(dlange,DLANGE) 51 #define DGECON_F77 F77_FUNC(dgecon,DGECON) 52 #define DGGLSE_F77 F77_FUNC(dgglse,DGGLSE) 53 #define DSYEV_F77 F77_FUNC(dsyev,DSYEV) 54 #define DSYTRF_F77 F77_FUNC(dsytrf,DSYTRF) 55 #define DSYTRI_F77 F77_FUNC(dsytri,DSYTRI) 56 #define DSYTRS_F77 F77_FUNC(dsytrs,DSYTRS) 57 #define DTRTRS_F77 F77_FUNC(dtrtrs,DTRTRS) 58 #define DSYCON_F77 F77_FUNC(dsycon,DSYCON) 59 // 60 #else 61 // Use the CMake generated fortran name mangling macros (eliminate warnings) 62 63 // BMA (2016): Changed to use Fortran 2003 ISO C bindings. 64 // The Fortran symbol will be lowercase with same name as if in C 65 //#define NKM_PIVOTCHOL_F77 SURF77_GLOBAL(nkm_pivotchol,NKM_PIVOTCHOL) 66 //#define NKM_BLOCKPIVOTCHOL_F77 SURF77_GLOBAL(nkm_blockpivotchol,NKM_BLOCKPIVOTCHOL) 67 #define NKM_PIVOTCHOL_F77 nkm_pivotchol 68 #define NKM_BLOCKPIVOTCHOL_F77 nkm_blockpivotchol 69 70 71 #define DGETRF_F77 SURF77_GLOBAL(dgetrf,DGETRF) 72 #define DGETRI_F77 SURF77_GLOBAL(dgetri,DGETRI) 73 #define DGEMV_F77 SURF77_GLOBAL(dgemv,DGEMV) 74 #define DGEMM_F77 SURF77_GLOBAL(dgemm,DGEMM) 75 #define DDOT_F77 SURF77_GLOBAL(ddot, DDOT) 76 #define DGELS_F77 SURF77_GLOBAL(dgels,DGELS) 77 #define DGESVD_F77 SURF77_GLOBAL(dgesvd,DGESVD) 78 79 #define DPOTRF_F77 SURF77_GLOBAL(dpotrf,DPOTRF) 80 #define DPOTRI_F77 SURF77_GLOBAL(dpotri,DPOTRI) 81 #define DPOTRS_F77 SURF77_GLOBAL(dpotrs,DPOTRS) 82 #define DPOCON_F77 SURF77_GLOBAL(dpocon,DPOCON) 83 #define DGETRS_F77 SURF77_GLOBAL(dgetrs,DGETRS) 84 #define DLANGE_F77 SURF77_GLOBAL(dlange,DLANGE) 85 #define DGECON_F77 SURF77_GLOBAL(dgecon,DGECON) 86 #define DGGLSE_F77 SURF77_GLOBAL(dgglse,DGGLSE) 87 #define DSYEV_F77 SURF77_GLOBAL(dsyev,DSYEV) 88 #define DSYTRF_F77 SURF77_GLOBAL(dsytrf,DSYTRF) 89 #define DSYTRI_F77 SURF77_GLOBAL(dsytri,DSYTRI) 90 #define DSYTRS_F77 SURF77_GLOBAL(dsytrs,DSYTRS) 91 #define DSYCON_F77 SURF77_GLOBAL(dsycon,DSYCON) 92 #define DTRTRS_F77 SURF77_GLOBAL(dtrtrs,DTRTRS) 93 // 94 #endif 95 96 /***************************************************************************/ 97 /**** Fortran to C name mangling ****/ 98 /***************************************************************************/ 99 100 extern "C" { // prevent C++ name mangling 101 102 //wrapper for KRD's implementation of the pivoting cholesky algorithm of: C. Lucas, "LAPACK-style Codes for LEvel2 and 3 Pivoted Cholesky Factorizations", Numerical Analysis Report No. 442, February 2004, from the Manchester Center for Computational Mathematics, I downloaded it from http://www.maths.manchester.ac.uk/~nareports/narep442.pdf, I think that Lucas's lev3pchol.f was ALWAYS defaulting to his lev2pchol.f, but this implementation was between 5% and 10% faster than Lucas's for the test problem (a 5500x5500 R matrix paviani 10D, 500 pts) 103 void NKM_PIVOTCHOL_F77(const char* uplo, const int* n, double* A, const int* lda, 104 int* piv, int* rank, const double* tol, int* info); 105 void NKM_BLOCKPIVOTCHOL_F77(const char* uplo, const int* n, double* A, 106 const int* lda, const int *blocksize, int* piv, 107 int* rank, double *dwork, const double* tol, 108 int* info); 109 110 111 /***************************************************************************/ 112 /**** BLAS Fortran to C name mangling ****/ 113 /***************************************************************************/ 114 115 // Vector-vector inner product 116 double DDOT_F77(const int* n, const double* x, const int* incx, 117 const double* y, const int* incy); 118 119 120 // Matrix-vector multiplication 121 void DGEMV_F77(char* trans, const int* m, const int* n, const double* alpha, 122 const double* A, const int* lda, const double* x, 123 const int* incx, const double* beta, double* y, const int* incy); 124 125 // Matrix-matrix multiplication 126 void DGEMM_F77(char* transa, char* transb, const int* m, const int* n, 127 const int* k, const double* alpha, const double* A, 128 const int* lda, const double* B, const int* ldb, 129 const double* beta, double* C, const int* ldc); 130 131 /***************************************************************************/ 132 /**** LAPACK Fortran to C name mangling ****/ 133 /***************************************************************************/ 134 135 // Perform Cholesky factorization 136 void DPOTRF_F77(const char* uplo, const int* n, double* AChol, const int* lda, int* info); 137 138 // Compute the inverse of a matrix expressed as an cholesky decomposition (i.e., call dpotrf on the matrix first) 139 void DPOTRI_F77(const char* uplo, const int* n, double* ACholInv, const int* lda, int* info); 140 141 // solve A*X=B for X, after A has been Cholesky factorized (i.e., call dptorf on the matrix first) 142 void DPOTRS_F77(const char* uplo, const int* n, const int* nRHS, 143 const double* AChol, 144 const int* ldAChol , double* RHS, 145 const int* ldRHS, int* info); 146 147 // function to compute the condition number of a matrix from the Cholesky factorization 148 void DPOCON_F77(const char* uplo, const int* n, const double* AChol, const int* lda, const double* anorm, double* rconda, double* work, int* iwork, int* info); 149 150 151 // Performs an L*D*L^T (or U*D*U^T, either can be used) factorization 152 void DSYTRF_F77(const char* uplo, const int* n, double* ALDLT, const int* lda, int* ipiv, double* work, const int* lwork, const int* info); 153 154 // Compute the inverse of a matrix expressed as an L*D*L^T (or U*D*U^T, either can be used) factorization (i.e. call dsytrf on the matrix first) 155 void DSYTRI_F77(const char* uplo, const int* n, double* ALDLTINV, const int* lda, const int* ipiv, double* work, int* info); 156 157 // solve A*X=B for X, after A has been L*D*L^T (or U*D*U^T, either can be used) factorized (i.e. call dsytrf on the matrix first) 158 void DSYTRS_F77(const char* uplo, const int* n, const int* nRHS, const double* ALDLT, const int* lda, const int* ipiv, double* RHS, const int* ldRHS, int* info); 159 160 // function to compute the condition number of a matrix from the L*D*L^T (or U*D*U^T, either can be used) factorization 161 void DSYCON_F77(const char* uplo, const int* n, const double* ALDLT, const int* lda, const int* ipiv, const double* anorm, double* rconda, double* work, int* iwork, int* info); 162 163 164 // Perform SVD factorization 165 void DGESVD_F77(const char* jobu, const char* jobvt, const int* m, const int* n, double* A, const int* lda, double* S, double* U, const int* ldu, double* VT, const int* ldvt, double* work, const int* lwork, int* info); 166 167 168 // Perform LU factorization 169 void DGETRF_F77(const int* m, const int* n, double* a, const int* lda, 170 int* ipiv, int* info); 171 172 // Compute the inverse of a matrix expressed as an LU decomposition (i.e., call dgetrf on the matrix first) 173 void DGETRI_F77(const int* n, double* a, const int* lda, const int* ipiv, 174 double* work, const int* lwork, int* info); 175 176 // solve A*X=B for X, where A={A || A^T} after A has been LU factorized (i.e., call dgetrf on the matrix first) 177 void DGETRS_F77(const char* transLU, const int* n, const int* nRHS, 178 const double* LU, const int* ldLU , const int* ipiv, 179 double* RHS, const int* ldRHS, int* info); 180 181 //function to compute the norm of a matrix A, choices are 182 //M max(abs(A(i,j))) this is not a consistent matrix norm, 183 //1 one norm of a matrix, maximum column sum, 184 //I infinity norm of matrix, maximum row sum,or 185 //F frobenius norm of a matrix, square root of sum of squares 186 double DLANGE_F77(char *whichnorm, int *M, int *N, const double *A, int *LDA, 187 double *work); 188 189 //function to compute the condition number of a matrix 190 void DGECON_F77(const char *whichnorm, const int *N, const double *ALU, 191 const int *LDA, const double *anorm, 192 double *rcond, double *work, int *iwork, int *info); 193 194 // Least-squares solution to linear system of equations 195 void DGELS_F77(const char* trans, const int* nrows, const int* ncols, 196 const int* nrhs, double* A, const int* lda, double* b, 197 const int* ldb, double* work, const int* lwork, int* info); 198 199 // Performs least-squares solve subject to equality constraints 200 void DGGLSE_F77(const int* m, const int* n, const int* p, double* A, 201 const int* lda, double* B, const int* ldb, double* c, 202 double* d, double* x, double* work, const int* lwork, 203 int* info); 204 205 // determines eigenvalues and (optionally) eigenvectors for a real symmetric matrix 206 void DSYEV_F77(const char* jobz, const char* uplo, const int* N, 207 double *A_EIGVECT, const int* lda, double* eigval, 208 double* work, const int* lwork, int* info); 209 210 // Computes the solution to AX=B when A is triangular 211 void DTRTRS_F77(const char* uplo, const char* trans, const char* diag, 212 const int* ncols, const int* nrhs, const double* A, 213 const int* lda, double* B, const int* ldb, int* info); 214 215 } // extern "C" (prevent C++ name mangling) 216 217 218 #endif // __SURFPACK_LAPACK_WRAPPERS_H__ 219 220