Line data Source code
1 : /** \file templateLapack.cpp
2 : * \brief Implementation of templatized wrappers for the Lapack library
3 : * \ingroup gen_math_files
4 : *
5 : */
6 :
7 : //***********************************************************************//
8 : // Copyright 2020 Jared R. Males (jaredmales@gmail.com)
9 : //
10 : // This file is part of mxlib.
11 : //
12 : // mxlib is free software: you can redistribute it and/or modify
13 : // it under the terms of the GNU General Public License as published by
14 : // the Free Software Foundation, either version 3 of the License, or
15 : // (at your option) any later version.
16 : //
17 : // mxlib is distributed in the hope that it will be useful,
18 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
19 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 : // GNU General Public License for more details.
21 : //
22 : // You should have received a copy of the GNU General Public License
23 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
24 : //***********************************************************************//
25 :
26 : #include "math/templateLapack.hpp"
27 :
28 : #ifndef MXLIB_MKL
29 : // Generic LAPACKE headers do not declare the provider's LAED9 Fortran symbols.
30 : extern "C"
31 : {
32 : void slaed9_( const MXLAPACK_INT *K,
33 : const MXLAPACK_INT *KSTART,
34 : const MXLAPACK_INT *KSTOP,
35 : const MXLAPACK_INT *N,
36 : float *D,
37 : float *Q,
38 : const MXLAPACK_INT *LDQ,
39 : const float *RHO,
40 : float *DLAMDA,
41 : float *W,
42 : float *S,
43 : const MXLAPACK_INT *LDS,
44 : MXLAPACK_INT *INFO );
45 :
46 : void dlaed9_( const MXLAPACK_INT *K,
47 : const MXLAPACK_INT *KSTART,
48 : const MXLAPACK_INT *KSTOP,
49 : const MXLAPACK_INT *N,
50 : double *D,
51 : double *Q,
52 : const MXLAPACK_INT *LDQ,
53 : const double *RHO,
54 : double *DLAMDA,
55 : double *W,
56 : double *S,
57 : const MXLAPACK_INT *LDS,
58 : MXLAPACK_INT *INFO );
59 : }
60 : #endif
61 :
62 : namespace mx
63 : {
64 : namespace math
65 : {
66 :
67 : template <>
68 11 : float lamch<float>( char CMACH )
69 : {
70 11 : return slamch_( &CMACH
71 : #ifdef LAPACK_FORTRAN_STRLEN_END
72 : ,
73 : 1
74 : #endif
75 11 : );
76 : }
77 :
78 : // Double specialization of lamch, a wrapper for Lapack DLAMCH
79 : template <>
80 171 : double lamch<double>( char CMACH )
81 : {
82 171 : return dlamch_( &CMACH
83 : #ifdef LAPACK_FORTRAN_STRLEN_END
84 : ,
85 : 1
86 : #endif
87 171 : );
88 : }
89 :
90 : template <>
91 0 : MXLAPACK_INT potrf<float>( char UPLO, MXLAPACK_INT N, float *A, MXLAPACK_INT LDA, MXLAPACK_INT &INFO )
92 : {
93 0 : spotrf_( &UPLO,
94 : &N,
95 : A,
96 : &LDA,
97 : &INFO
98 : #ifdef LAPACK_FORTRAN_STRLEN_END
99 : ,
100 : 1
101 : #endif
102 : );
103 :
104 0 : return INFO;
105 : }
106 :
107 : template <>
108 0 : MXLAPACK_INT potrf<double>( char UPLO, MXLAPACK_INT N, double *A, MXLAPACK_INT LDA, MXLAPACK_INT &INFO )
109 : {
110 0 : dpotrf_( &UPLO,
111 : &N,
112 : A,
113 : &LDA,
114 : &INFO
115 : #ifdef LAPACK_FORTRAN_STRLEN_END
116 : ,
117 : 1
118 : #endif
119 : );
120 :
121 0 : return INFO;
122 : }
123 :
124 : template <>
125 : MXLAPACK_INT
126 0 : potrf<std::complex<float>>( char UPLO, MXLAPACK_INT N, std::complex<float> *A, MXLAPACK_INT LDA, MXLAPACK_INT &INFO )
127 : {
128 0 : cpotrf_( &UPLO,
129 : &N,
130 : (float _Complex *)A,
131 : &LDA,
132 : &INFO
133 : #ifdef LAPACK_FORTRAN_STRLEN_END
134 : ,
135 : 1
136 : #endif
137 : );
138 :
139 0 : return INFO;
140 : }
141 :
142 : template <>
143 : MXLAPACK_INT
144 0 : potrf<std::complex<double>>( char UPLO, MXLAPACK_INT N, std::complex<double> *A, MXLAPACK_INT LDA, MXLAPACK_INT &INFO )
145 : {
146 0 : zpotrf_( &UPLO,
147 : &N,
148 : (double _Complex *)A,
149 : &LDA,
150 : &INFO
151 : #ifdef LAPACK_FORTRAN_STRLEN_END
152 : ,
153 : 1
154 : #endif
155 : );
156 :
157 0 : return INFO;
158 : }
159 :
160 : template <>
161 0 : MXLAPACK_INT sytrd<float>( char UPLO,
162 : MXLAPACK_INT N,
163 : float *A,
164 : MXLAPACK_INT LDA,
165 : float *D,
166 : float *E,
167 : float *TAU,
168 : float *WORK,
169 : MXLAPACK_INT LWORK,
170 : MXLAPACK_INT INFO )
171 : {
172 :
173 0 : ssytrd_( &UPLO,
174 : &N,
175 : A,
176 : &LDA,
177 : D,
178 : E,
179 : TAU,
180 : WORK,
181 : &LWORK,
182 : &INFO
183 : #ifdef LAPACK_FORTRAN_STRLEN_END
184 : ,
185 : 1
186 : #endif
187 : );
188 :
189 0 : return INFO;
190 : }
191 :
192 : template <>
193 0 : MXLAPACK_INT sytrd<double>( char UPLO,
194 : MXLAPACK_INT N,
195 : double *A,
196 : MXLAPACK_INT LDA,
197 : double *D,
198 : double *E,
199 : double *TAU,
200 : double *WORK,
201 : MXLAPACK_INT LWORK,
202 : MXLAPACK_INT INFO )
203 : {
204 :
205 0 : dsytrd_( &UPLO,
206 : &N,
207 : A,
208 : &LDA,
209 : D,
210 : E,
211 : TAU,
212 : WORK,
213 : &LWORK,
214 : &INFO
215 : #ifdef LAPACK_FORTRAN_STRLEN_END
216 : ,
217 : 1
218 : #endif
219 : );
220 :
221 0 : return INFO;
222 : }
223 :
224 : template <>
225 2 : MXLAPACK_INT laed9<float>( float *D,
226 : float *Q,
227 : float *S,
228 : MXLAPACK_INT K,
229 : MXLAPACK_INT KSTART,
230 : MXLAPACK_INT KSTOP,
231 : MXLAPACK_INT N,
232 : MXLAPACK_INT LDQ,
233 : float RHO,
234 : float *DLAMDA,
235 : float *W,
236 : MXLAPACK_INT LDS )
237 : {
238 : MXLAPACK_INT INFO;
239 :
240 2 : slaed9_( &K, &KSTART, &KSTOP, &N, D, Q, &LDQ, &RHO, DLAMDA, W, S, &LDS, &INFO );
241 :
242 2 : return INFO;
243 : }
244 :
245 : template <>
246 11 : MXLAPACK_INT laed9<double>( double *D,
247 : double *Q,
248 : double *S,
249 : MXLAPACK_INT K,
250 : MXLAPACK_INT KSTART,
251 : MXLAPACK_INT KSTOP,
252 : MXLAPACK_INT N,
253 : MXLAPACK_INT LDQ,
254 : double RHO,
255 : double *DLAMDA,
256 : double *W,
257 : MXLAPACK_INT LDS )
258 : {
259 : MXLAPACK_INT INFO;
260 :
261 11 : dlaed9_( &K, &KSTART, &KSTOP, &N, D, Q, &LDQ, &RHO, DLAMDA, W, S, &LDS, &INFO );
262 :
263 11 : return INFO;
264 : }
265 :
266 : // Float specialization of syevr, a wrapper for Lapack SSYEVR
267 : template <>
268 10 : MXLAPACK_INT syevr<float>( char JOBZ,
269 : char RANGE,
270 : char UPLO,
271 : MXLAPACK_INT N,
272 : float *A,
273 : MXLAPACK_INT LDA,
274 : float VL,
275 : float VU,
276 : MXLAPACK_INT IL,
277 : MXLAPACK_INT IU,
278 : float ABSTOL,
279 : MXLAPACK_INT *M,
280 : float *W,
281 : float *Z,
282 : MXLAPACK_INT LDZ,
283 : MXLAPACK_INT *ISUPPZ,
284 : float *WORK,
285 : MXLAPACK_INT LWORK,
286 : MXLAPACK_INT *IWORK,
287 : MXLAPACK_INT LIWORK )
288 : {
289 :
290 : MXLAPACK_INT INFO;
291 :
292 10 : ssyevr_( &JOBZ,
293 : &RANGE,
294 : &UPLO,
295 : &N,
296 : A,
297 : &LDA,
298 : &VL,
299 : &VU,
300 : &IL,
301 : &IU,
302 : &ABSTOL,
303 : M,
304 : W,
305 : Z,
306 : &LDZ,
307 : ISUPPZ,
308 : WORK,
309 : &LWORK,
310 : IWORK,
311 : &LIWORK,
312 : &INFO
313 : #ifdef LAPACK_FORTRAN_STRLEN_END
314 : ,
315 : 1,
316 : 1,
317 : 1
318 : #endif
319 : );
320 :
321 10 : return INFO;
322 : }
323 :
324 : // Double specialization of syevr, a wrapper for Lapack DSYEVR
325 : template <>
326 137 : MXLAPACK_INT syevr<double>( char JOBZ,
327 : char RANGE,
328 : char UPLO,
329 : MXLAPACK_INT N,
330 : double *A,
331 : MXLAPACK_INT LDA,
332 : double VL,
333 : double VU,
334 : MXLAPACK_INT IL,
335 : MXLAPACK_INT IU,
336 : double ABSTOL,
337 : MXLAPACK_INT *M,
338 : double *W,
339 : double *Z,
340 : MXLAPACK_INT LDZ,
341 : MXLAPACK_INT *ISUPPZ,
342 : double *WORK,
343 : MXLAPACK_INT LWORK,
344 : MXLAPACK_INT *IWORK,
345 : MXLAPACK_INT LIWORK )
346 : {
347 :
348 : MXLAPACK_INT INFO;
349 :
350 137 : dsyevr_( &JOBZ,
351 : &RANGE,
352 : &UPLO,
353 : &N,
354 : A,
355 : &LDA,
356 : &VL,
357 : &VU,
358 : &IL,
359 : &IU,
360 : &ABSTOL,
361 : M,
362 : W,
363 : Z,
364 : &LDZ,
365 : ISUPPZ,
366 : WORK,
367 : &LWORK,
368 : IWORK,
369 : &LIWORK,
370 : &INFO
371 : #ifdef LAPACK_FORTRAN_STRLEN_END
372 : ,
373 : 1,
374 : 1,
375 : 1
376 : #endif
377 : );
378 :
379 137 : return INFO;
380 : }
381 :
382 : // float specialization of gesvd
383 : template <>
384 5 : MXLAPACK_INT gesvd<float>( char JOBU,
385 : char JOBVT,
386 : MXLAPACK_INT M,
387 : MXLAPACK_INT N,
388 : float *A,
389 : MXLAPACK_INT LDA,
390 : float *S,
391 : float *U,
392 : MXLAPACK_INT LDU,
393 : float *VT,
394 : MXLAPACK_INT LDVT,
395 : float *WORK,
396 : MXLAPACK_INT LWORK )
397 : {
398 : MXLAPACK_INT INFO;
399 :
400 5 : sgesvd_( &JOBU,
401 : &JOBVT,
402 : &M,
403 : &N,
404 : A,
405 : &LDA,
406 : S,
407 : U,
408 : &LDU,
409 : VT,
410 : &LDVT,
411 : WORK,
412 : &LWORK,
413 : &INFO
414 : #ifdef LAPACK_FORTRAN_STRLEN_END
415 : ,
416 : 1,
417 : 1
418 : #endif
419 : );
420 :
421 5 : return INFO;
422 : }
423 :
424 : // double specialization of gesvd
425 : template <>
426 55 : MXLAPACK_INT gesvd<double>( char JOBU,
427 : char JOBVT,
428 : MXLAPACK_INT M,
429 : MXLAPACK_INT N,
430 : double *A,
431 : MXLAPACK_INT LDA,
432 : double *S,
433 : double *U,
434 : MXLAPACK_INT LDU,
435 : double *VT,
436 : MXLAPACK_INT LDVT,
437 : double *WORK,
438 : MXLAPACK_INT LWORK )
439 : {
440 : MXLAPACK_INT INFO;
441 :
442 55 : dgesvd_( &JOBU,
443 : &JOBVT,
444 : &M,
445 : &N,
446 : A,
447 : &LDA,
448 : S,
449 : U,
450 : &LDU,
451 : VT,
452 : &LDVT,
453 : WORK,
454 : &LWORK,
455 : &INFO
456 : #ifdef LAPACK_FORTRAN_STRLEN_END
457 : ,
458 : 1,
459 : 1
460 : #endif
461 : );
462 :
463 55 : return INFO;
464 : }
465 :
466 : // float specialization of gesdd
467 : template <>
468 0 : MXLAPACK_INT gesdd<float>( char JOBZ,
469 : MXLAPACK_INT M,
470 : MXLAPACK_INT N,
471 : float *A,
472 : MXLAPACK_INT LDA,
473 : float *S,
474 : float *U,
475 : MXLAPACK_INT LDU,
476 : float *VT,
477 : MXLAPACK_INT LDVT,
478 : float *WORK,
479 : MXLAPACK_INT LWORK,
480 : MXLAPACK_INT *IWORK )
481 : {
482 : MXLAPACK_INT INFO;
483 :
484 0 : sgesdd_( &JOBZ,
485 : &M,
486 : &N,
487 : A,
488 : &LDA,
489 : S,
490 : U,
491 : &LDU,
492 : VT,
493 : &LDVT,
494 : WORK,
495 : &LWORK,
496 : IWORK,
497 : &INFO
498 : #ifdef LAPACK_FORTRAN_STRLEN_END
499 : ,
500 : 1
501 : #endif
502 : );
503 :
504 0 : return INFO;
505 : }
506 :
507 : // double specialization of gesdd
508 : template <>
509 0 : MXLAPACK_INT gesdd<double>( char JOBZ,
510 : MXLAPACK_INT M,
511 : MXLAPACK_INT N,
512 : double *A,
513 : MXLAPACK_INT LDA,
514 : double *S,
515 : double *U,
516 : MXLAPACK_INT LDU,
517 : double *VT,
518 : MXLAPACK_INT LDVT,
519 : double *WORK,
520 : MXLAPACK_INT LWORK,
521 : MXLAPACK_INT *IWORK )
522 : {
523 :
524 : MXLAPACK_INT INFO;
525 :
526 0 : dgesdd_( &JOBZ,
527 : &M,
528 : &N,
529 : A,
530 : &LDA,
531 : S,
532 : U,
533 : &LDU,
534 : VT,
535 : &LDVT,
536 : WORK,
537 : &LWORK,
538 : IWORK,
539 : &INFO
540 : #ifdef LAPACK_FORTRAN_STRLEN_END
541 : ,
542 : 1
543 : #endif
544 : );
545 :
546 0 : return INFO;
547 : }
548 :
549 : } // namespace math
550 : } // namespace mx
|
