Line data Source code
1 : !-----------------------------------------------------------------------
2 : ! $Id$
3 : !===============================================================================
4 : module lapack_interfaces
5 :
6 : ! This module acts as an interface to Lapack routines. It's main purpose
7 : ! is to make interfaces available to clubb that can handle both
8 : ! single and doulbe precision. This may be compiled along with Lapack
9 : ! source code, or along with a linked Lapack library such as MKL.
10 :
11 : implicit none
12 :
13 : public :: lapack_gbsv, lapack_gbsvx, &
14 : lapack_gtsv, lapack_gtsvx, &
15 : lapack_isnan, lapack_potrf, &
16 : lapack_poequ, lapack_laqsy, &
17 : lapack_syev, lapack_trmv
18 :
19 : private :: &
20 : dgbsv_wrap, sgbsv_wrap, &
21 : dgbsvx_wrap, sgbsvx_wrap, &
22 : dgtsv_wrap, sgtsv_wrap, &
23 : dgtsvx_wrap, sgtsvx_wrap, &
24 : disnan_wrap, sisnan_wrap, &
25 : dpotrf_wrap, spotrf_wrap, &
26 : dpoequ_wrap, spoequ_wrap, &
27 : dlaqsy_wrap, slaqsy_wrap, &
28 : dsyev_wrap, ssyev_wrap, &
29 : dtrmv_wrap, strmv_wrap
30 :
31 : ! Interface for Lapack general band solver, single or double precision
32 : interface lapack_gbsv
33 : module procedure dgbsv_wrap
34 : module procedure sgbsv_wrap
35 : end interface
36 :
37 : ! Interface for Lapack general band solver, expert version, single or double precision
38 : interface lapack_gbsvx
39 : module procedure dgbsvx_wrap
40 : module procedure sgbsvx_wrap
41 : end interface
42 :
43 : ! Interface for Lapack tridiagonal matrix solver, single or double precision
44 : interface lapack_gtsv
45 : module procedure dgtsv_wrap
46 : module procedure sgtsv_wrap
47 : end interface
48 :
49 : ! Interface for Lapack tridiagonal matrix solver, expert version, single or double precision
50 : interface lapack_gtsvx
51 : module procedure dgtsvx_wrap
52 : module procedure sgtsvx_wrap
53 : end interface
54 :
55 : ! Interface for Lapack nan check, single or double precision
56 : interface lapack_isnan
57 : module procedure disnan_wrap
58 : module procedure sisnan_wrap
59 : end interface
60 :
61 : ! Interface for Lapack's Cholesky factorization of a real symmetric positive definite
62 : ! matrix, single or double precision
63 : interface lapack_potrf
64 : module procedure dpotrf_wrap
65 : module procedure spotrf_wrap
66 : end interface
67 :
68 : ! Interface for Lapack routine to compute row and column scalings intended to
69 : ! equilibriate a symmetric positive definite matrix, single or doulbe precision
70 : interface lapack_poequ
71 : module procedure dpoequ_wrap
72 : module procedure spoequ_wrap
73 : end interface
74 :
75 : ! Interface for Lapack routine to equilibriate a symmetric matrix, single or double precision
76 : interface lapack_laqsy
77 : module procedure dlaqsy_wrap
78 : module procedure slaqsy_wrap
79 : end interface
80 :
81 : ! Interface for Lapack routine to compute all eigenvalues and, optionally, eigenvectors
82 : ! of a real symmetric matrix, single or double precision
83 : interface lapack_syev
84 : module procedure dsyev_wrap
85 : module procedure ssyev_wrap
86 : end interface
87 :
88 : ! Interface for Lapack routines to performe one of the following matrix-vector operations
89 : ! x := A*x, or x := A**T*x,
90 : ! where A is an upper or lower triangular matrix, single or double precision
91 : interface lapack_trmv
92 : module procedure dtrmv_wrap
93 : module procedure strmv_wrap
94 : end interface
95 :
96 :
97 : private ! Set Default Scope
98 :
99 : contains
100 :
101 : ! ==================== General Band Solver Wrappers ====================
102 :
103 : ! Double precision wrapper
104 280519200 : subroutine dgbsv_wrap( n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info )
105 :
106 : implicit none
107 :
108 : external :: dgbsv
109 :
110 : integer info, kl, ku, ldab, ldb, n, nrhs
111 : integer ipiv( * )
112 : double precision ab( ldab, * ), b( ldb, * )
113 :
114 280519200 : call dgbsv( n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info )
115 :
116 280519200 : end subroutine dgbsv_wrap
117 :
118 : ! Single precision wrapper
119 0 : subroutine sgbsv_wrap( n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info )
120 :
121 : implicit none
122 :
123 : external :: sgbsv
124 :
125 : integer info, kl, ku, ldab, ldb, n, nrhs
126 : integer ipiv( * )
127 : real ab( ldab, * ), b( ldb, * )
128 :
129 0 : call sgbsv( n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info )
130 :
131 0 : end subroutine sgbsv_wrap
132 :
133 :
134 : ! ==================== Band Solver Expert Wrappers ====================
135 :
136 : ! Double precision wrapper
137 0 : subroutine dgbsvx_wrap( fact, trans, n, kl, ku, nrhs, ab, ldab, afb, &
138 0 : ldafb, ipiv, equed, r, c, b, ldb, x, ldx, &
139 : rcond, ferr, berr, work, iwork, info )
140 : implicit none
141 :
142 : external :: dgbsvx
143 :
144 : character equed, fact, trans
145 : integer info, kl, ku, ldab, ldafb, ldb, ldx, n, nrhs
146 : double precision rcond
147 : integer ipiv( * ), iwork( * )
148 : double precision ab( ldab, * ), afb( ldafb, * ), b( ldb, * ), &
149 : berr( * ), c( * ), ferr( * ), r( * ), &
150 : work( * ), x( ldx, * )
151 :
152 : call dgbsvx( fact, trans, n, kl, ku, nrhs, ab, ldab, afb, &
153 : ldafb, ipiv, equed, r, c, b, ldb, x, ldx, &
154 0 : rcond, ferr, berr, work, iwork, info )
155 :
156 0 : end subroutine dgbsvx_wrap
157 :
158 : ! Single precision wrapper
159 0 : subroutine sgbsvx_wrap( fact, trans, n, kl, ku, nrhs, ab, ldab, afb, &
160 0 : ldafb, ipiv, equed, r, c, b, ldb, x, ldx, &
161 : rcond, ferr, berr, work, iwork, info )
162 : implicit none
163 :
164 : external :: sgbsvx
165 :
166 : character equed, fact, trans
167 : integer info, kl, ku, ldab, ldafb, ldb, ldx, n, nrhs
168 : real rcond
169 : integer ipiv( * ), iwork( * )
170 : real ab( ldab, * ), afb( ldafb, * ), b( ldb, * ), &
171 : berr( * ), c( * ), ferr( * ), r( * ), &
172 : work( * ), x( ldx, * )
173 :
174 : call sgbsvx( fact, trans, n, kl, ku, nrhs, ab, ldab, afb, &
175 : ldafb, ipiv, equed, r, c, b, ldb, x, ldx, &
176 0 : rcond, ferr, berr, work, iwork, info )
177 :
178 0 : end subroutine sgbsvx_wrap
179 :
180 :
181 : ! ==================== Tridiagonal Solver Wrappers ====================
182 :
183 : ! Double precision wrapper
184 707664624 : subroutine dgtsv_wrap( n, nrhs, dl, d, du, b, ldb, info )
185 :
186 : implicit none
187 :
188 : external :: dgtsv
189 :
190 : integer info, ldb, n, nrhs
191 : double precision b( ldb, * ), d( * ), dl( * ), du( * )
192 :
193 707664624 : call dgtsv( n, nrhs, dl, d, du, b, ldb, info )
194 :
195 707664624 : end subroutine dgtsv_wrap
196 :
197 : ! Single precision wrapper
198 0 : subroutine sgtsv_wrap( n, nrhs, dl, d, du, b, ldb, info )
199 :
200 : implicit none
201 :
202 : external :: sgtsv
203 :
204 : integer info, ldb, n, nrhs
205 : real b( ldb, * ), d( * ), dl( * ), du( * )
206 :
207 0 : call sgtsv( n, nrhs, dl, d, du, b, ldb, info )
208 :
209 0 : end subroutine sgtsv_wrap
210 :
211 :
212 : ! ==================== Tridiagonal Solver Expert Wrappers ====================
213 :
214 : ! Double precision wrapper
215 0 : subroutine dgtsvx_wrap( fact, trans, n, nrhs, dl, d, du, dlf, df, duf, &
216 0 : du2, ipiv, b, ldb, x, ldx, rcond, ferr, berr, &
217 : work, iwork, info )
218 : implicit none
219 :
220 : external :: dgtsvx
221 :
222 : character fact, trans
223 : integer info, ldb, ldx, n, nrhs
224 : double precision rcond
225 : integer ipiv( * ), iwork( * )
226 : double precision b( ldb, * ), berr( * ), d( * ), df( * ), &
227 : dl( * ), dlf( * ), du( * ), du2( * ), duf( * ), &
228 : ferr( * ), work( * ), x( ldx, * )
229 :
230 : call dgtsvx( fact, trans, n, nrhs, dl, d, du, dlf, df, duf, &
231 : du2, ipiv, b, ldb, x, ldx, rcond, ferr, berr, &
232 0 : work, iwork, info )
233 :
234 0 : end subroutine dgtsvx_wrap
235 :
236 : ! Single precision wrapper
237 0 : subroutine sgtsvx_wrap( fact, trans, n, nrhs, dl, d, du, dlf, df, duf, &
238 0 : du2, ipiv, b, ldb, x, ldx, rcond, ferr, berr, &
239 : work, iwork, info )
240 : implicit none
241 :
242 : external :: sgtsvx
243 :
244 : character fact, trans
245 : integer info, ldb, ldx, n, nrhs
246 : real rcond
247 : integer ipiv( * ), iwork( * )
248 : real b( ldb, * ), berr( * ), d( * ), df( * ), &
249 : dl( * ), dlf( * ), du( * ), du2( * ), duf( * ), &
250 : ferr( * ), work( * ), x( ldx, * )
251 :
252 : call sgtsvx( fact, trans, n, nrhs, dl, d, du, dlf, df, duf, &
253 : du2, ipiv, b, ldb, x, ldx, rcond, ferr, berr, &
254 0 : work, iwork, info )
255 :
256 0 : end subroutine sgtsvx_wrap
257 :
258 :
259 : ! ==================== NaN Check Wrappers ====================
260 :
261 : ! Double precision wrapper
262 : !-----------------------------------------------------------------------
263 707664624 : logical function disnan_wrap( ndim, nrhs, variable )
264 :
265 : ! Description:
266 : ! Check for NaN values in a variable using the LAPACK subroutines
267 :
268 : ! References:
269 : ! <http://www.netlib.org/lapack/single/sisnan.f>
270 : ! <http://www.netlib.org/lapack/double/disnan.f>
271 : !-----------------------------------------------------------------------
272 :
273 : implicit none
274 :
275 : #ifdef NO_LAPACK_ISNAN /* Used for older LAPACK libraries that don't have sisnan/disnan */
276 :
277 : integer, intent(in) :: &
278 : ndim, & ! Size of variable
279 : nrhs ! Number of right hand sides
280 :
281 : double precision, dimension(ndim,nrhs), intent(in) :: &
282 : variable ! Variable to check
283 :
284 >21837*10^7 : disnan_wrap = any( variable(:,1:nrhs) /= variable(:,1:nrhs) )
285 : #else
286 :
287 : logical, external :: &
288 : disnan ! Procedure
289 :
290 : integer, intent(in) :: &
291 : ndim, & ! Size of variable
292 : nrhs ! Number of right hand sides
293 :
294 : double precision, dimension(ndim,nrhs), intent(in) :: &
295 : variable ! Variable to check
296 :
297 : integer :: k, j
298 :
299 : ! ---- Begin Code ----
300 :
301 : disnan_wrap = .false.
302 :
303 : do k = 1, ndim
304 : do j = 1, nrhs
305 :
306 : ! Lapack NaN check function, sisnan for single precision or disnan for double precision
307 : disnan_wrap = disnan( variable(k,j) )
308 :
309 : if ( disnan_wrap ) exit
310 : end do
311 : if ( disnan_wrap ) exit
312 : end do
313 :
314 : #endif /* NO_LAPACK_ISNAN */
315 :
316 : return
317 : end function disnan_wrap
318 :
319 : ! Single precision wrapper
320 : !-----------------------------------------------------------------------
321 0 : logical function sisnan_wrap( ndim, nrhs, variable )
322 :
323 : ! Description:
324 : ! Check for NaN values in a variable using the LAPACK subroutines
325 :
326 : ! References:
327 : ! <http://www.netlib.org/lapack/single/sisnan.f>
328 : ! <http://www.netlib.org/lapack/double/disnan.f>
329 : !-----------------------------------------------------------------------
330 :
331 : implicit none
332 :
333 : #ifdef NO_LAPACK_ISNAN /* Used for older LAPACK libraries that don't have sisnan/disnan */
334 :
335 : integer, intent(in) :: &
336 : ndim, & ! Size of variable
337 : nrhs ! Number of right hand sides
338 :
339 : real, dimension(ndim,nrhs), intent(in) :: &
340 : variable ! Variable to check
341 :
342 0 : sisnan_wrap = any( variable(:,1:nrhs) /= variable(:,1:nrhs) )
343 : #else
344 :
345 : logical, external :: &
346 : sisnan ! Procedure
347 :
348 : integer, intent(in) :: &
349 : ndim, & ! Size of variable
350 : nrhs ! Number of right hand sides
351 :
352 : real, dimension(ndim,nrhs), intent(in) :: &
353 : variable ! Variable to check
354 :
355 : integer :: k, j
356 :
357 : ! ---- Begin Code ----
358 :
359 : sisnan_wrap = .false.
360 :
361 : do k = 1, ndim
362 : do j = 1, nrhs
363 :
364 : ! Lapack NaN check function, sisnan for single precision or disnan for double precision
365 : sisnan_wrap = sisnan( variable(k,j) )
366 :
367 : if ( sisnan_wrap ) exit
368 : end do
369 : if ( sisnan_wrap ) exit
370 : end do
371 :
372 : #endif /* NO_LAPACK_ISNAN */
373 :
374 : return
375 : end function sisnan_wrap
376 :
377 :
378 : ! ==================== Cholesky Factorization Wrappers ====================
379 :
380 : ! Double precision wrapper
381 0 : subroutine dpotrf_wrap( uplo, n, a, lda, info )
382 :
383 : implicit none
384 :
385 : external :: dpotrf
386 :
387 : character uplo
388 : integer info, lda, n
389 : double precision a( lda, * )
390 :
391 0 : call dpotrf( uplo, n, a, lda, info )
392 :
393 0 : end subroutine dpotrf_wrap
394 :
395 : ! Single precision wrapper
396 0 : subroutine spotrf_wrap( uplo, n, a, lda, info )
397 :
398 : implicit none
399 :
400 : external :: spotrf
401 :
402 : character uplo
403 : integer info, lda, n
404 : real a( lda, * )
405 :
406 0 : call spotrf( uplo, n, a, lda, info )
407 :
408 0 : end subroutine spotrf_wrap
409 :
410 :
411 : ! ==================== Equilibrium Scaling Calculation Wrappers ====================
412 :
413 : ! Double precision wrapper
414 0 : subroutine dpoequ_wrap( n, a, lda, s, scond, amax, info )
415 :
416 : implicit none
417 :
418 : external :: dpoequ
419 :
420 : integer info, lda, n
421 : double precision amax, scond
422 : double precision a( lda, * ), s( * )
423 :
424 0 : call dpoequ( n, a, lda, s, scond, amax, info )
425 :
426 0 : end subroutine dpoequ_wrap
427 :
428 : ! Single precision wrapper
429 0 : subroutine spoequ_wrap( n, a, lda, s, scond, amax, info )
430 :
431 : implicit none
432 :
433 : external :: spoequ
434 :
435 : integer info, lda, n
436 : real amax, scond
437 : real a( lda, * ), s( * )
438 :
439 0 : call spoequ( n, a, lda, s, scond, amax, info )
440 :
441 0 : end subroutine spoequ_wrap
442 :
443 :
444 : ! ==================== Matrix Equilibriator Wrappers ====================
445 :
446 : ! Double precision wrapper
447 0 : subroutine dlaqsy_wrap( uplo, n, a, lda, s, scond, amax, equed )
448 :
449 : implicit none
450 :
451 : external :: dlaqsy
452 :
453 : character equed, uplo
454 : integer lda, n
455 : double precision amax, scond
456 : double precision a( lda, * ), s( * )
457 :
458 0 : call dlaqsy( uplo, n, a, lda, s, scond, amax, equed )
459 :
460 0 : end subroutine dlaqsy_wrap
461 :
462 : ! Single precision wrapper
463 0 : subroutine slaqsy_wrap( uplo, n, a, lda, s, scond, amax, equed )
464 :
465 : implicit none
466 :
467 : external :: slaqsy
468 :
469 : character equed, uplo
470 : integer lda, n
471 : real amax, scond
472 : real a( lda, * ), s( * )
473 :
474 0 : call slaqsy( uplo, n, a, lda, s, scond, amax, equed )
475 :
476 0 : end subroutine slaqsy_wrap
477 :
478 :
479 : ! ==================== Eigenvalue/vector Calculation Wrappers ====================
480 :
481 : ! Double precision wrapper
482 0 : subroutine dsyev_wrap( jobz, uplo, n, a, lda, w, work, lwork, info )
483 :
484 : implicit none
485 :
486 : external :: dsyev
487 :
488 : character jobz, uplo
489 : integer info, lda, lwork, n
490 : double precision a( lda, * ), w( * ), work( * )
491 :
492 0 : call dsyev( jobz, uplo, n, a, lda, w, work, lwork, info )
493 :
494 0 : end subroutine dsyev_wrap
495 :
496 : ! Single precision wrapper
497 0 : subroutine ssyev_wrap( jobz, uplo, n, a, lda, w, work, lwork, info )
498 :
499 : implicit none
500 :
501 : external :: ssyev
502 :
503 : character jobz, uplo
504 : integer info, lda, lwork, n
505 : real a( lda, * ), w( * ), work( * )
506 :
507 0 : call ssyev( jobz, uplo, n, a, lda, w, work, lwork, info )
508 :
509 0 : end subroutine ssyev_wrap
510 :
511 :
512 : ! ==================== Matrix Operations Wrappers ====================
513 :
514 : ! Double precision wrapper
515 0 : subroutine dtrmv_wrap( uplo, trans, diag, n, a, lda, x, incx)
516 :
517 : implicit none
518 :
519 : external :: dtrmv
520 :
521 : integer incx,lda,n
522 : character diag,trans,uplo
523 : double precision a(lda,*),x(*)
524 :
525 0 : call dtrmv( uplo, trans, diag, n, a, lda, x, incx)
526 :
527 0 : end subroutine dtrmv_wrap
528 :
529 : ! Single precision wrapper
530 0 : subroutine strmv_wrap( uplo, trans, diag, n, a, lda, x, incx)
531 :
532 : implicit none
533 :
534 : external :: strmv
535 :
536 : integer incx,lda,n
537 : character diag,trans,uplo
538 : real a(lda,*),x(*)
539 :
540 0 : call strmv( uplo, trans, diag, n, a, lda, x, incx)
541 :
542 0 : end subroutine strmv_wrap
543 :
544 : end module lapack_interfaces
|