1 /*  $Id: ncbi_lapack.c 503534 2016-06-06 14:44:06Z ucko $
2 * ===========================================================================
3 *
4 *                            PUBLIC DOMAIN NOTICE
5 *               National Center for Biotechnology Information
6 *
7 *  This software/database is a "United States Government Work" under the
8 *  terms of the United States Copyright Act.  It was written as part of
9 *  the author's official duties as a United States Government employee and
10 *  thus cannot be copyrighted.  This software/database is freely available
11 *  to the public for use. The National Library of Medicine and the U.S.
12 *  Government have not placed any restriction on its use or reproduction.
13 *
14 *  Although all reasonable efforts have been taken to ensure the accuracy
15 *  and reliability of the software and data, the NLM and the U.S.
16 *  Government do not and cannot warrant the performance or results that
17 *  may be obtained by using this software or data. The NLM and the U.S.
18 *  Government disclaim all warranties, express or implied, including
19 *  warranties of performance, merchantability or fitness for any particular
20 *  purpose.
21 *
22 *  Please cite the author in any work or product based on this material.
23 *
24 * ===========================================================================
25 *
26 * Author:  Aaron Ucko
27 *
28 * File Description:
29 *   Wrappers for LAPACK routines of interest.
30 *
31 * ===========================================================================
32 */
33 
34 #include <misc/lapackwrapp/ncbi_lapack.h>
35 #include <stdlib.h>
36 #include <string.h>
37 
38 #include <ncbiconf.h>
39 
40 #ifdef HAVE_LAPACKE_H
41 #  define HAVE_LAPACK_CONFIG_H 1
42 #  include <lapacke.h>
43 #elif defined(HAVE_LAPACKE_LAPACKE_H)
44 #  define HAVE_LAPACK_CONFIG_H 1
45 #  include <lapacke/lapacke.h>
46 #elif defined(HAVE_ACCELERATE_ACCELERATE_H)
47 #  include <Accelerate/Accelerate.h>
48 #elif defined(HAVE___CLPK_INTEGER)
49 #  include <clapack.h>
50 #else
51 typedef int                     TLapackInt;
52 typedef int                     TLapackLogical;
53 typedef struct { float  r, i; } TLapackComplexFloat;
54 typedef struct { double r, i; } TLapackComplexDouble;
55 #  ifdef HAVE_CLAPACK_H
56 #    define complex       TLapackComplexFloat
57 #    define doublecomplex TLapackComplexDouble
58 #    define integer       TLapackInt
59 #    define logical       TLapackLogical
60 #    define VOID          void
61 #    include <clapack.h>
62 #    undef complex
63 #    undef doublecomplex
64 #    undef integer
65 #    undef logical
66 #    undef VOID
67 #  endif
68 #endif
69 
70 #if defined(LAPACK_malloc)
71 typedef lapack_int            TLapackInt;
72 typedef lapack_logical        TLapackLogical;
73 typedef lapack_complex_float  TLapackComplexFloat;
74 typedef lapack_complex_double TLapackComplexDouble;
75 #elif defined(HAVE___CLPK_INTEGER)
76 typedef __CLPK_integer       TLapackInt;
77 typedef __CLPK_logical       TLapackLogical;
78 typedef __CLPK_complex       TLapackComplexFloat;
79 typedef __CLPK_doublecomplex TLapackComplexDouble;
80 #endif
81 
82 /* Fall back as necessary to supplying our own prototypes to work around
83  * https://bugzilla.redhat.com/show_bug.cgi?id=1165538 or the like. */
84 #if !defined(LAPACK_COL_MAJOR)  &&  !defined(__CLAPACK_H)
85 TLapackInt dsyev_(char* jobz, char* uplo, TLapackInt* n,
86                   double* a, TLapackInt* lda, double* w,
87                   double* work, TLapackInt* lwork, TLapackInt* info);
88 TLapackInt dgels_(char* trans, TLapackInt* m, TLapackInt* n,
89                   TLapackInt* nrhs, double* a, TLapackInt* lda,
90                   double* b, TLapackInt* ldb,
91                   double* work, TLapackInt* lwork, TLapackInt* info);
92 TLapackInt dgesvd_(char* jobu, char* jobvt, TLapackInt* m, TLapackInt* n,
93                    double* a, TLapackInt* lda, double* s,
94                    double* u, TLapackInt* ldu, double* vt, TLapackInt* ldvt,
95                    double* work, TLapackInt* lwork, TLapackInt* info);
96 #endif
97 
s_Min(unsigned int x,unsigned int y)98 static unsigned int s_Min(unsigned int x, unsigned int y)
99 {
100     return x < y ? x : y;
101 }
102 
NCBI_SymmetricEigens(enum EEigensWanted jobz_in,enum EMatrixTriangle uplo_in,unsigned int n_in,double * a,unsigned int lda_in,double * w)103 int NCBI_SymmetricEigens(enum EEigensWanted jobz_in,
104                          enum EMatrixTriangle uplo_in, unsigned int n_in,
105                          double* a, unsigned int lda_in, double* w)
106 {
107     double     optimal_work_size;
108     char       jobz = jobz_in, uplo = uplo_in;
109     TLapackInt n = n_in, lda = lda_in, lwork = -1, info = 0;
110     dsyev_(&jobz, &uplo, &n, a, &lda, w, &optimal_work_size, &lwork, &info);
111     if (info == 0) {
112         double* work;
113         lwork = optimal_work_size;
114         work  = malloc(lwork * sizeof(double));
115         if (work == NULL) {
116             return kNcbiLapackMemoryError;
117         }
118         dsyev_(&jobz, &uplo, &n, a, &lda, w, work, &lwork, &info);
119         free(work);
120     }
121     return info;
122 }
123 
NCBI_LinearSolution(enum EMaybeTransposed trans_in,unsigned int r,unsigned int c,unsigned int nrhs_in,double * a,unsigned int lda_in,double * b,unsigned int ldb_in)124 int NCBI_LinearSolution(enum EMaybeTransposed trans_in,
125                         unsigned int r, unsigned int c, unsigned int nrhs_in,
126                         double* a, unsigned int lda_in,
127                         double* b, unsigned int ldb_in)
128 {
129     double     ows;
130     char       trans = trans_in;
131     TLapackInt m = r, n = c, nrhs = nrhs_in, lda = lda_in, ldb = ldb_in,
132                lwork = -1, info = 0;
133     dgels_(&trans, &m, &n, &nrhs, a, &lda, b, &ldb, &ows, &lwork, &info);
134     if (info == 0) {
135         double* work;
136         lwork = ows;
137         work  = malloc(lwork * sizeof(double));
138         if (work == NULL) {
139             return kNcbiLapackMemoryError;
140         }
141         dgels_(&trans, &m, &n, &nrhs, a, &lda, b, &ldb, work, &lwork, &info);
142         free(work);
143     }
144     return info;
145 }
146 
NCBI_SingularValueDecomposition(enum ESingularVectorDetailsWanted jobu_in,enum ESingularVectorDetailsWanted jobvt_in,unsigned int r,unsigned int c,double * a,unsigned int lda_in,double * s,double * u,unsigned int ldu_in,double * vt,unsigned int ldvt_in,double * superb)147 int NCBI_SingularValueDecomposition(enum ESingularVectorDetailsWanted jobu_in,
148                                     enum ESingularVectorDetailsWanted jobvt_in,
149                                     unsigned int r, unsigned int c,
150                                     double* a, unsigned int lda_in, double* s,
151                                     double* u, unsigned int ldu_in,
152                                     double* vt, unsigned int ldvt_in,
153                                     double* superb)
154 {
155     double     optimal_work_size;
156     char       jobu = jobu_in, jobvt = jobvt_in;
157     TLapackInt m = r, n = c, lda = lda_in, ldu = ldu_in, ldvt = ldvt_in,
158                lwork = -1, info = 0;
159     dgesvd_(&jobu, &jobvt, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt,
160             &optimal_work_size, &lwork, &info);
161     if (info == 0) {
162         double* work;
163         lwork = optimal_work_size;
164         work  = malloc(lwork * sizeof(double));
165         if (work == NULL) {
166             return kNcbiLapackMemoryError;
167         }
168         dgesvd_(&jobu, &jobvt, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt,
169                 work, &lwork, &info);
170         if (superb != NULL /* && info > 0 */) {
171             memcpy(superb, work + 1, (s_Min(r, c) - 1) * sizeof(double));
172         }
173         free(work);
174     }
175     return info;
176 }
177