ROMS
Loading...
Searching...
No Matches
lapack_mod.F
Go to the documentation of this file.
1#include "cppdefs.h"
2 MODULE lapack_mod
3!
4!git $Id$
5!================================================== Hernan G. Arango ===
6! Copyright (c) 2002-2025 The ROMS Group !
7! Licensed under a MIT/X style license !
8! See License_ROMS.md !
9!=======================================================================
10! !
11! This module includes modernized versions of the selected routines !
12! from thw Linear Algebra Package (LAPACK) library, which are used !
13! in ROMS 4D-Var algorithms. !
14! !
15! Adapted from LAPACK library version 2.0 !
16! !
17! NOTES: !
18! !
19! - The LAPACK library was written originally in Fortran-77 and has !
20! been extensively tested and used in various compilers and !
21! applications. We are modernizing a few of the functions used by !
22! ROMS at the request of NOAA and NCEP to remove the undesirable !
23! GOTOs statements with the modern capabilities of the Fortran !
24! standard (1995 and 2003). !
25! !
26! - We need to declare the non-scalar arguments to routines as !
27! assumed-size array, A(*), instead of assumed-shape array, A(:). !
28! It is relatively common in linear algebra to pass sub-matrices !
29! as an argument. However, this feature was declared obsolescent !
30! in Fortran95. The colon we need an explicit interface. Currently,!
31! if fails to compile if colon are used intead of asterisk, even !
32! if it is inside this module. !
33! !
34!======================================================================!
35!
36 USE mod_kinds
37!
38! implicit none
39!
40 PUBLIC :: dsteqr
41!
42 PRIVATE :: dlae2
43 PRIVATE :: dlaev2
44 PRIVATE :: dlamc1
45 PRIVATE :: dlamc2
46 PRIVATE :: dlamc3
47 PRIVATE :: dlamc4
48 PRIVATE :: dlamc5
49 PRIVATE :: dlamch
50 PRIVATE :: dlanst
51 PRIVATE :: dlapy2
52 PRIVATE :: dlartg
53 PRIVATE :: dlascl
54 PRIVATE :: dlaset
55 PRIVATE :: dlasr
56 PRIVATE :: dlasrt
57 PRIVATE :: dlassq
58 PRIVATE :: lsame
59 PRIVATE :: xerbla
60!
61!-----------------------------------------------------------------------
62 CONTAINS
63!-----------------------------------------------------------------------
64!
65 SUBROUTINE dsteqr ( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
66!
67!=======================================================================
68! !
69! DSTEQR computes all eigenvalues and, optionally, eigenvectors of a !
70! symmetric tridiagonal matrix using the implicit QL or QR method. !
71! The eigenvectors of a full or band symmetric matrix can also be !
72! found if DSYTRD or DSPTRD or DSBTRD has been used to reduce this !
73! matrix to tridiagonal form. !
74! !
75! Arguments: !
76! !
77! COMPZ (input) CHARACTER*1 !
78! = 'N': Compute eigenvalues only. !
79! = 'V': Compute eigenvalues and eigenvectors of the !
80! original symmetric matrix. On entry, Z must contain !
81! the orthogonal matrix used to reduce the original !
82! matrix to tridiagonal form. !
83! = 'I': Compute eigenvalues and eigenvectors of the !
84! tridiagonal matrix. Z is initialized to the !
85! identity matrix. !
86! !
87! N (input) INTEGER !
88! The order of the matrix. N >= 0. !
89! !
90! D (input/output) DOUBLE PRECISION array, dimension (N) !
91! On entry, the diagonal elements of the tridiagonal matrix. !
92! On exit, if INFO = 0, the eigenvalues in ascending order. !
93! !
94! E (input/output) DOUBLE PRECISION array, dimension (N-1) !
95! On entry, the (n-1) subdiagonal elements of the tridiagonal !
96! matrix. !
97! On exit, E has been destroyed. !
98! !
99! Z (input/output) DOUBLE PRECISION array, dimension (LDZ, N) !
100! On entry, if COMPZ = 'V', then Z contains the orthogonal !
101! matrix used in the reduction to tridiagonal form. !
102! On exit, if INFO = 0, then if COMPZ = 'V', Z contains the !
103! orthonormal eigenvectors of the original symmetric matrix, !
104! and if COMPZ = 'I', Z contains the orthonormal eigenvectors !
105! of the symmetric tridiagonal matrix. !
106! If COMPZ = 'N', then Z is not referenced. !
107! !
108! LDZ (input) INTEGER !
109! The leading dimension of the array Z. LDZ >= 1, and if !
110! eigenvectors are desired, then LDZ >= max(1,N). !
111! !
112! WORK (workspace) DOUBLE PRECISION array, dimension !
113! (max(1,2*N-2)). !
114! If COMPZ = 'N', then WORK is not referenced. !
115! !
116! INFO (output) INTEGER !
117! = 0: successful exit !
118! < 0: if INFO = -i, the i-th argument had an illegal value !
119! > 0: the algorithm has failed to find all the eigenvalues !
120! in a total of 30*N iterations; if INFO = i, then i !
121! elements of E have not converged to zero; on exit, D !
122! and E contain the elements of a symmetric tridiagonal !
123! matrix which is orthogonally similar to the original !
124! matrix. !
125! !
126! -- LAPACK routine (version 2.0) -- !
127! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., !
128! Courant Institute, Argonne National Lab, and Rice University !
129! September 30, 1994 !
130! !
131!=======================================================================
132!
133! Imported variable declarations.
134!
135 integer, intent(in) :: ldz, n
136 integer, intent(out) :: info
137!
138 real (dp), intent(inout) :: work(*)
139 real (dp), intent(inout) :: d(*), e(*), z(ldz,*)
140!
141 character (len=1), intent(in) :: compz
142!
143! Local variable declarations.
144!
145 integer, parameter :: maxit = 30
146!
147 real (dp), parameter :: zero = 0.0_dp
148 real (dp), parameter :: one = 1.0_dp
149 real (dp), parameter :: two = 2.0_dp
150 real (dp), parameter :: three = 3.0_dp
151!
152 integer :: i, icompz, ii, iscale, j, jtot, k, l, l1, lend, lendm1
153 integer :: lendp1, lendsv, lm1, lsv, m, mm, mm1, nm1, nmaxit
154!
155 real (dp) :: anorm, b, c, eps, eps2, f, g, p, r, rt1, rt2
156 real (dp) :: s, safmax, safmin, ssfmax, ssfmin, tst
157!
158!-----------------------------------------------------------------------
159! Executable Statements.
160!-----------------------------------------------------------------------
161!
162! Test the input parameters.
163!
164 info = 0
165!
166 IF (lsame(compz, 'N')) THEN
167 icompz = 0
168 ELSE IF (lsame(compz, 'V')) THEN
169 icompz = 1
170 ELSE IF (lsame(compz, 'I' )) THEN
171 icompz = 2
172 ELSE
173 icompz = -1
174 END IF
175!
176 IF (icompz.LT.0) THEN
177 info = -1
178 ELSE IF( n.LT.0 ) THEN
179 info = -2
180 ELSE IF ((ldz.LT.1).OR.(icompz.GT.0 .AND. ldz.LT.max(1,n))) THEN
181 info = -6
182 END IF
183!
184 IF (info.NE.0) THEN
185 CALL xerbla( 'DSTEQR', -info )
186 RETURN
187 END IF
188!
189! Quick return if possible
190!
191 IF (n.EQ.0) RETURN
192!
193 IF ( n.EQ.1 ) THEN
194 IF (icompz.EQ.2) z(1,1) = one
195 RETURN
196 END IF
197!
198! Determine the unit roundoff and over/underflow thresholds.
199!
200 eps = dlamch('E')
201 eps2 = eps**2
202 safmin = dlamch('S')
203 safmax = one / safmin
204 ssfmax = sqrt(safmax) / three
205 ssfmin = sqrt(safmin) / eps2
206!
207! Compute the eigenvalues and eigenvectors of the tridiagonal
208! matrix.
209!
210 IF (icompz .EQ. 2) THEN
211 CALL dlaset ('Full', n, n, zero, one, z, ldz)
212 END IF
213!
214 nmaxit = n*maxit
215 jtot = 0
216!
217! Determine where the matrix splits and choose QL or QR iteration
218! for each block, according to whether top or bottom diagonal
219! element is smaller.
220!
221 l1 = 1
222 nm1 = n - 1
223!
224 10 CONTINUE
225 IF (l1 .GT. n) GO TO 160
226 IF (l1 .GT. 1) e(l1-1) = zero
227 IF (l1 .LE. nm1) THEN
228 DO m = l1, nm1
229 tst = abs(e(m))
230 IF (tst .EQ. zero) GO TO 30
231 IF (tst .LE. (sqrt(abs(d(m)))*sqrt(abs(d(m+1))))*eps) THEN
232 e(m) = zero
233 GO TO 30
234 END IF
235 END DO
236 END IF
237 m = n
238!
239 30 CONTINUE
240 l = l1
241 lsv = l
242 lend = m
243 lendsv = lend
244 l1 = m + 1
245 IF (lend .EQ. l) GO TO 10
246!
247! Scale submatrix in rows and columns L to LEND.
248!
249 anorm = dlanst('I', lend-l+1, d(l), e(l))
250 iscale = 0
251 IF (anorm .EQ. zero) GO TO 10
252 IF (anorm .GT. ssfmax) THEN
253 iscale = 1
254 CALL dlascl ('G', 0, 0, anorm, ssfmax, lend-l+1, 1, d(l), n, &
255 & info )
256 CALL dlascl ('G', 0, 0, anorm, ssfmax, lend-l, 1, e(l), n, &
257 & info )
258 ELSE IF (anorm .LT. ssfmin) THEN
259 iscale = 2
260 CALL dlascl ('G', 0, 0, anorm, ssfmin, lend-l+1, 1, d(l), n, &
261 & info )
262 CALL dlascl ('G', 0, 0, anorm, ssfmin, lend-l, 1, e( l ), n, &
263 & info )
264 END IF
265!
266! Choose between QL and QR iteration
267!
268 IF (abs(d(lend)) .LT. abs(d(l))) THEN
269 lend = lsv
270 l = lendsv
271 END IF
272!
273 IF (lend .GT. l) THEN
274!
275! QL Iteration.
276!
277! Look for small subdiagonal element.
278!
279 40 CONTINUE
280 IF (l .NE. lend) THEN
281 lendm1 = lend - 1
282 DO m = l, lendm1
283 tst = abs(e(m))**2
284 IF (tst .LE. (eps2*abs(d(m)))*abs(d(m+1))+safmin) GO TO 60
285 END DO
286 END IF
287!
288 m = lend
289!
290 60 CONTINUE
291 IF (m .LT. lend) e(m) = zero
292 p = d(l)
293 IF (m .EQ. l) GO TO 80
294!
295! If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
296! to compute its eigensystem.
297!
298 IF (m .EQ. l+1) THEN
299 IF (icompz .GT. 0) THEN
300 CALL dlaev2 ( d(l), e(l), d(l+1), rt1, rt2, c, s)
301 work(l) = c
302 work(n-1+l) = s
303 CALL dlasr ('R', 'V', 'B', n, 2, work(l), &
304 & work(n-1+l), z(1, l), ldz)
305 ELSE
306 CALL dlae2 (d(l), e(l), d(l+1), rt1, rt2)
307 END IF
308 d(l) = rt1
309 d(l+1) = rt2
310 e(l) = zero
311 l = l + 2
312 IF (l .LE. lend) GO TO 40
313 GO TO 140
314 END IF
315!
316 IF (jtot .EQ. nmaxit) GO TO 140
317 jtot = jtot + 1
318!
319! Form shift.
320!
321 g = (d(l+1)-p) / (two*e(l))
322 r = dlapy2(g, one)
323 g = d(m) - p + (e(l) / (g + sign(r, g)))
324!
325 s = one
326 c = one
327 p = zero
328!
329! Inner loop
330!
331 mm1 = m - 1
332 DO i = mm1, l, -1
333 f = s*e(i)
334 b = c*e(i)
335 CALL dlartg (g, f, c, s, r)
336 IF (i .NE. m-1) e(i+1) = r
337 g = d(i+1) - p
338 r = (d(i)-g)*s + two*c*b
339 p = s*r
340 d(i+1) = g+p
341 g = c*r - b
342!
343! If eigenvectors are desired, then save rotations.
344!
345 IF (icompz.GT.0) THEN
346 work(i) = c
347 work(n-1+i) = -s
348 END IF
349 END DO
350!
351! If eigenvectors are desired, then apply saved rotations.
352!
353 IF (icompz.GT.0) THEN
354 mm = m - l + 1
355 CALL dlasr ('R', 'V', 'B', n, mm, work(l), work(n-1+l), &
356 & z(1,l), ldz)
357 END IF
358!
359 d(l) = d( l ) - p
360 e(l) = g
361 GO TO 40
362!
363! Eigenvalue found.
364!
365 80 CONTINUE
366 d(l) = p
367 !
368 l = l + 1
369 IF (l .LE. lend) GO TO 40
370 GO TO 140
371!
372 ELSE
373!
374! QR Iteration
375!
376! Look for small superdiagonal element.
377!
378 90 CONTINUE
379 IF (l .NE. lend) THEN
380 lendp1 = lend + 1
381 DO m = l, lendp1, -1
382 tst = abs(e( m-1))**2
383 IF (tst .LE. (eps2*abs(d(m)))*abs(d(m-1))+safmin) GO TO 110
384 END DO
385 END IF
386!
387 m = lend
388!
389 110 CONTINUE
390 IF (m .GT. lend) e(m-1) = zero
391 p = d(l)
392 IF (m .EQ. l) GO TO 130
393!
394! If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
395! to compute its eigensystem.
396!
397 IF (m .EQ. l-1) THEN
398 IF (icompz .GT. 0) THEN
399 CALL dlaev2 (d(l-1), e(l-1), d(l), rt1, rt2, c, s)
400 work(m) = c
401 work(n-1+m) = s
402 CALL dlasr ('R', 'V', 'F', n, 2, work(m), &
403 & work(n-1+m), z(1,l-1), ldz)
404 ELSE
405 CALL dlae2 (d(l-1), e(l-1), d(l), rt1, rt2)
406 END IF
407 d(l-1) = rt1
408 d(l) = rt2
409 e(l-1) = zero
410 l = l - 2
411 IF (l .GE. lend) GO TO 90
412 GO TO 140
413 END IF
414!
415 IF (jtot .EQ. nmaxit) GO TO 140
416 jtot = jtot + 1
417!
418! Form shift.
419!
420 g = (d(l-1)-p) / (two*e(l-1))
421 r = dlapy2(g, one)
422 g = d(m) - p + (e(l-1) / (g+sign(r, g)))
423!
424 s = one
425 c = one
426 p = zero
427!
428! Inner loop
429!
430 lm1 = l - 1
431 DO i = m, lm1
432 f = s*e(i)
433 b = c*e(i)
434 CALL dlartg (g, f, c, s, r)
435 IF (i .NE. m) e(i-1) = r
436 g = d(i) - p
437 r = (d(i+1)-g)*s + two*c*b
438 p = s*r
439 d(i) = g + p
440 g = c*r - b
441!
442! If eigenvectors are desired, then save rotations.
443!
444 IF (icompz.GT.0) THEN
445 work(i) = c
446 work(n-1+i) = s
447 END IF
448 END DO
449!
450! If eigenvectors are desired, then apply saved rotations.
451!
452 IF (icompz.GT.0) THEN
453 mm = l - m + 1
454 CALL dlasr ('R', 'V', 'F', n, mm, work(m), work(n-1+m), &
455 & z(1,m), ldz)
456 END IF
457!
458 d(l) = d(l) - p
459 e(lm1) = g
460 GO TO 90
461!
462! Eigenvalue found.
463!
464 130 CONTINUE
465 d(l) = p
466!
467 l = l - 1
468 IF (l .GE. lend) GO TO 90
469 GO TO 140
470!
471 END IF
472!
473! Undo scaling if necessary.
474!
475 140 CONTINUE
476 IF (iscale .EQ. 1) THEN
477 CALL dlascl ('G', 0, 0, ssfmax, anorm, lendsv-lsv+1, 1, &
478 & d(lsv), n, info)
479 CALL dlascl ('G', 0, 0, ssfmax, anorm, lendsv-lsv, 1, e(lsv), &
480 & n, info)
481 ELSE IF (iscale .EQ. 2) THEN
482 CALL dlascl ('G', 0, 0, ssfmin, anorm, lendsv-lsv+1, 1, &
483 & d(lsv), n, info)
484 CALL dlascl ('G', 0, 0, ssfmin, anorm, lendsv-lsv, 1, e(lsv), &
485 & n, info)
486 END IF
487!
488! Check for no convergence to an eigenvalue after a total
489! of N*MAXIT iterations.
490!
491 IF (jtot .LT. nmaxit) GO TO 10
492 DO i = 1, n - 1
493 IF (e(i) .NE. zero) info = info + 1
494 END DO
495 GO TO 190
496!
497! Order eigenvalues and eigenvectors.
498!
499 160 CONTINUE
500 IF (icompz .EQ. 0) THEN
501!
502! Use Quick Sort.
503!
504 CALL dlasrt ('I', n, d, info)
505!
506 ELSE
507!
508! Use Selection Sort to minimize swaps of eigenvectors.
509!
510 DO ii = 2, n
511 i = ii - 1
512 k = i
513 p = d( i )
514 DO j = ii, n
515 IF (d(j) .LT. p) THEN
516 k = j
517 p = d(j)
518 END IF
519 END DO
520 IF (k .NE. i) THEN
521 d(k) = d(i)
522 d(i) = p
523 CALL dswap (n, z(1,i), 1, z(1,k), 1)
524 END IF
525 END DO
526 END IF
527!
528 190 CONTINUE
529 RETURN
530!
531 END SUBROUTINE dsteqr
532!
533 SUBROUTINE dlae2 (A, B, C, RT1, RT2)
534!
535!=======================================================================
536! !
537! DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix: !
538! !
539! [ A B ] !
540! [ B C ]. !
541! !
542! On return, RT1 is the eigenvalue of larger absolute value, and RT2 !
543! is the eigenvalue of smaller absolute value. !
544! !
545! Arguments: !
546! !
547! A (input) DOUBLE PRECISION !
548! The (1,1) element of the 2-by-2 matrix. !
549! !
550! B (input) DOUBLE PRECISION !
551! The (1,2) and (2,1) elements of the 2-by-2 matrix. !
552! !
553! C (input) DOUBLE PRECISION !
554! The (2,2) element of the 2-by-2 matrix. !
555! !
556! RT1 (output) DOUBLE PRECISION !
557! The eigenvalue of larger absolute value. !
558! !
559! RT2 (output) DOUBLE PRECISION !
560! The eigenvalue of smaller absolute value. !
561! !
562! Further Details: !
563! !
564! RT1 is accurate to a few ulps barring over/underflow. !
565! !
566! RT2 may be inaccurate if there is massive cancellation in the !
567! determinant A*C-B*B; higher precision or correctly rounded or !
568! correctly truncated arithmetic would be needed to compute RT2 !
569! accurately in all cases. !
570! !
571! Overflow is possible only if RT1 is within a factor of 5 of !
572! overflow. !
573! !
574! Underflow is harmless if the input data is 0 or exceeds !
575! underflow_threshold / macheps. !
576! !
577! -- LAPACK auxiliary routine (version 2.0) -- !
578! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., !
579! Courant Institute, Argonne National Lab, and Rice University !
580! October 31, 1992 !
581! !
582!=======================================================================
583!
584! Imported variable declarations.
585!
586 real (dp), intent(in ) :: a, b, c
587 real (dp), intent(out) :: rt1, rt2
588!
589! Local variable declarations.
590!
591 real (dp), parameter :: one = 1.0_dp
592 real (dp), parameter :: two = 2.0_dp
593 real (dp), parameter :: zero = 0.0_dp
594 real (dp), parameter :: half = 0.5_dp
595!
596 real (dp) :: ab, acmn, acmx, adf, df, rt, sm, tb
597!
598!-----------------------------------------------------------------------
599! Executable Statements.
600!-----------------------------------------------------------------------
601!
602! Compute the eigenvalues.
603!
604 sm = a + c
605 df = a - c
606 adf = abs(df)
607 tb = b + b
608 ab = abs(tb)
609!
610 IF (abs(a) .GT. abs(c)) THEN
611 acmx = a
612 acmn = c
613 ELSE
614 acmx = c
615 acmn = a
616 END IF
617!
618 IF (adf .GT. ab) THEN
619 rt = adf*sqrt(one+(ab/adf )**2)
620 ELSE IF (adf .LT. ab) THEN
621 rt = ab*sqrt(one+(adf/ab )**2 )
622 ELSE
623 rt = ab*sqrt(two) ! Includes case AB=ADF=0
624 END IF
625!
626 IF (sm .LT. zero) THEN
627 rt1 = half*(sm-rt)
628!
629! Order of execution important.
630! To get fully accurate smaller eigenvalue,
631! next line needs to be executed in higher precision.
632!
633 rt2 = (acmx / rt1)*acmn - (b / rt1)*b
634 ELSE IF (sm .GT. zero ) THEN
635 rt1 = half*(sm + rt)
636!
637! Order of execution important.
638! To get fully accurate smaller eigenvalue,
639! next line needs to be executed in higher precision.
640!
641 rt2 = (acmx / rt1)*acmn - (b / rt1)*b
642 ELSE
643 rt1 = half*rt ! Includes case RT1 = RT2 = 0
644 rt2 = -half*rt
645 END IF
646!
647 RETURN
648 END SUBROUTINE dlae2
649!
650 SUBROUTINE dlaev2 (A, B, C, RT1, RT2, CS1, SN1)
651!
652!=======================================================================
653! !
654! DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric !
655! matrix: !
656! !
657! [ A B ] !
658! [ B C ]. !
659! !
660! On return, RT1 is the eigenvalue of larger absolute value, RT2 is !
661! the eigenvalue of smaller absolute value, and (CS1,SN1) is the unit !
662! right eigenvector for RT1, giving the decomposition: !
663! !
664! [ CS1 SN1 ] [ A B ] [ CS1 -SN1 ] = [ RT1 0 ] !
665! [-SN1 CS1 ] [ B C ] [ SN1 CS1 ] [ 0 RT2 ]. !
666! !
667! Arguments: !
668! !
669! A (input) DOUBLE PRECISION !
670! The (1,1) element of the 2-by-2 matrix. !
671! !
672! B (input) DOUBLE PRECISION !
673! The (1,2) element and the conjugate of the (2,1) element of !
674! the 2-by-2 matrix. !
675! !
676! C (input) DOUBLE PRECISION !
677! The (2,2) element of the 2-by-2 matrix. !
678! !
679! RT1 (output) DOUBLE PRECISION !
680! The eigenvalue of larger absolute value. !
681! !
682! RT2 (output) DOUBLE PRECISION !
683! The eigenvalue of smaller absolute value. !
684! !
685! CS1 (output) DOUBLE PRECISION !
686! SN1 (output) DOUBLE PRECISION !
687! The vector (CS1, SN1) is a unit right eigenvector for RT1. !
688! !
689! Further Details: !
690! !
691! RT1 is accurate to a few ulps barring over/underflow. !
692! !
693! RT2 may be inaccurate if there is massive cancellation in the !
694! determinant A*C-B*B; higher precision or correctly rounded or !
695! correctly truncated arithmetic would be needed to compute RT2 !
696! accurately in all cases. !
697! !
698! CS1 and SN1 are accurate to a few ulps barring over/underflow. !
699! !
700! Overflow is possible only if RT1 is within a factor of 5 of !
701! overflow. !
702! !
703! Underflow is harmless if the input data is 0 or exceeds !
704! underflow_threshold / macheps. !
705! !
706! -- LAPACK auxiliary routine (version 2.0) -- !
707! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., !
708! Courant Institute, Argonne National Lab, and Rice University !
709! October 31, 1992 !
710! !
711!=======================================================================
712!
713! Imported variable declarations.
714!
715 real (dp), intent(in ) :: a, b, c
716 real (dp), intent(out) :: cs1, rt1, rt2, sn1
717!
718! Local variable declarations.
719!
720 real (dp), parameter :: one = 1.0_dp
721 real (dp), parameter :: two = 2.0_dp
722 real (dp), parameter :: zero = 0.0_dp
723 real (dp), parameter :: half = 0.5_dp
724!
725 integer :: sgn1, sgn2
726 real (dp) :: ab, acmn, acmx, acs, adf, cs, ct, df, rt, sm, tb, tn
727!
728!-----------------------------------------------------------------------
729! Executable Statements.
730!-----------------------------------------------------------------------
731!
732! Compute the eigenvalues.
733!
734 sm = a + c
735 df = a - c
736 adf = abs(df)
737 tb = b + b
738 ab = abs(tb)
739!
740 IF (abs(a) .GT. abs(c)) THEN
741 acmx = a
742 acmn = c
743 ELSE
744 acmx = c
745 acmn = a
746 END IF
747!
748 IF (adf .GT. ab) THEN
749 rt = adf*sqrt(one+(ab/adf)**2)
750 ELSE IF (adf .LT. ab) THEN
751 rt = ab*sqrt(one+(adf/ab)**2 )
752 ELSE
753 rt = ab*sqrt(two) ! Includes case AB=ADF=0
754 END IF
755 IF (sm.LT.zero) THEN
756 rt1 = half*(sm - rt)
757 sgn1 = -1
758!
759! Order of execution important.
760! To get fully accurate smaller eigenvalue,
761! next line needs to be executed in higher precision.
762!
763 rt2 = (acmx/rt1)*acmn - (b/rt1)*b
764 ELSE IF (sm .GT. zero) THEN
765 rt1 = half*(sm + rt)
766 sgn1 = 1
767!
768! Order of execution important.
769! To get fully accurate smaller eigenvalue,
770! next line needs to be executed in higher precision.
771!
772 rt2 = (acmx/rt1)*acmn - (b/rt1)*b
773 ELSE
774 rt1 = half*rt ! Includes case RT1 = RT2 = 0
775 rt2 = -half*rt
776 sgn1 = 1
777 END IF
778!
779! Compute the eigenvector.
780!
781 IF (df .GE. zero) THEN
782 cs = df + rt
783 sgn2 = 1
784 ELSE
785 cs = df - rt
786 sgn2 = -1
787 END IF
788 acs = abs(cs)
789 IF (acs .GT. ab) THEN
790 ct = -tb / cs
791 sn1 = one / sqrt(one+ct*ct)
792 cs1 = ct*sn1
793 ELSE
794 IF (ab .EQ. zero) THEN
795 cs1 = one
796 sn1 = zero
797 ELSE
798 tn = -cs / tb
799 cs1 = one / sqrt(one+tn*tn)
800 sn1 = tn*cs1
801 END IF
802 END IF
803 IF (sgn1 .EQ. sgn2) THEN
804 tn = cs1
805 cs1 = -sn1
806 sn1 = tn
807 END IF
808!
809 RETURN
810 END SUBROUTINE dlaev2
811!
812 SUBROUTINE dlamc1 (BETA, T, RND, IEEE1)
813!
814!=======================================================================
815! !
816! DLAMC1 determines the machine parameters given by BETA, T, RND, and !
817! IEEE1. !
818! !
819! Arguments: !
820! !
821! BETA (output) INTEGER !
822! The base of the machine. !
823! !
824! T (output) INTEGER !
825! The number of ( BETA ) digits in the mantissa. !
826! !
827! RND (output) LOGICAL !
828! Specifies whether proper rounding ( RND = .TRUE. ) or !
829! chopping ( RND = .FALSE. ) occurs in addition. This may !
830! not be a reliable guide to the way in which the machine !
831! performs its arithmetic. !
832! !
833! IEEE1 (output) LOGICAL !
834! Specifies whether rounding appears to be done in the IEEE !
835! 'round to nearest' style. !
836! !
837! Further Details: !
838! !
839! The routine is based on the routine ENVRON by Malcolm and !
840! incorporates suggestions by Gentleman and Marovich. See !
841! !
842! Malcolm M. A. (1972) Algorithms to reveal properties of !
843! floating-point arithmetic. Comms. of the ACM, 15, 949-951. !
844! !
845! Gentleman W. M. and Marovich S. B. (1974) More on algorithms !
846! that reveal properties of floating point arithmetic units. !
847! Comms. of the ACM, 17, 276-277. !
848! !
849! -- LAPACK auxiliary routine (version 2.0) -- !
850! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., !
851! Courant Institute, Argonne National Lab, and Rice University !
852! October 31, 1992 !
853! !
854!=======================================================================
855!
856! Imported variable declarations.
857!
858 logical, intent(out) :: ieee1, rnd
859!
860 integer, intent(out) :: beta, t
861!
862! Local variable declarations.
863!
864 logical, save :: first = .true.
865 logical, save :: lieee1, lrnd
866!
867 integer, save :: lbeta, lt
868!
869 real (dp) :: a, b, c, f, one, qtr, savec, t1, t2
870!
871!-----------------------------------------------------------------------
872! Executable Statements.
873!-----------------------------------------------------------------------
874!
875 first_pass : IF (first) THEN
876 first = .false.
877 one = 1
878!
879! LBETA, LIEEE1, LT and LRND are the local values of BETA,
880! IEEE1, T and RND.
881!
882! Throughout this routine we use the function DLAMC3 to ensure
883! that relevant values are stored and not held in registers, or
884! are not affected by optimizers.
885!
886! Compute a = 2.0**m with the smallest positive integer m such that
887!
888! fl( a + 1.0 ) = a.
889!
890 a = 1
891 c = 1
892!
893 DO WHILE (c .EQ. one)
894 a = 2*a
895 c = dlamc3( a, one )
896 c = dlamc3( c, -a )
897 END DO
898!
899! Now compute b = 2.0**m with the smallest positive integer m
900! such that
901!
902! fl( a + b ) .gt. a.
903!
904 b = 1
905 c = dlamc3(a, b)
906!
907 DO WHILE (c .EQ. a)
908 b = 2*b
909 c = dlamc3(a, b)
910 END DO
911!
912! Now compute the base. a and c are neighbouring floating point
913! numbers in the interval ( beta**t, beta**( t + 1 ) ) and so
914! their difference is beta. Adding 0.25 to c is to ensure that it
915! is truncated to beta and not ( beta - 1 ).
916!
917 qtr = one / 4
918 savec = c
919 c = dlamc3(c, -a)
920 lbeta = c + qtr
921!
922! Now determine whether rounding or chopping occurs, by adding a
923! bit less than beta/2 and a bit more than beta/2 to a.
924!
925 b = lbeta
926 f = dlamc3(b/2, -b/100)
927 c = dlamc3(f, a)
928 IF (c .EQ. a) THEN
929 lrnd = .true.
930 ELSE
931 lrnd = .false.
932 END IF
933 f = dlamc3(b/2, b/100)
934 c = dlamc3(f, a)
935 IF ((lrnd) .AND. (c.EQ.a)) lrnd = .false.
936!
937! Try and decide whether rounding is done in the IEEE round to
938! nearest style. B/2 is half a unit in the last place of the two
939! numbers A and SAVEC. Furthermore, A is even, i.e. has last bit
940! zero, and SAVEC is odd. Thus adding B/2 to A should not change
941! A, but adding B/2 to SAVEC should change SAVEC.
942!
943 t1 = dlamc3(b/2, a)
944 t2 = dlamc3(b/2, savec)
945 lieee1 = (t1.EQ.a) .AND. (t2.GT.savec) .AND. lrnd
946!
947! Now find the mantissa, t. It should be the integer part of
948! log to the base beta of a, however it is safer to determine t
949! by powering. So we find t as the smallest positive integer for
950! which
951!
952! fl( beta**t + 1.0 ) = 1.0.
953!
954 lt = 0
955 a = 1
956 c = 1
957!
958 DO WHILE (c .EQ. one)
959 lt = lt + 1
960 a = a*lbeta
961 c = dlamc3(a, one)
962 c = dlamc3(c, -a)
963 END DO
964!
965 END IF first_pass
966!
967 beta = lbeta
968 t = lt
969 rnd = lrnd
970 ieee1 = lieee1
971!
972 RETURN
973 END SUBROUTINE dlamc1
974!
975 SUBROUTINE dlamc2 (BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX)
976!
977!=======================================================================
978! !
979! DLAMC2 determines the machine parameters specified in its argument !
980! list. !
981! !
982! Arguments: !
983! !
984! BETA (output) INTEGER !
985! The base of the machine. !
986! !
987! T (output) INTEGER !
988! The number of ( BETA ) digits in the mantissa. !
989! !
990! RND (output) LOGICAL !
991! Specifies whether proper rounding ( RND = .TRUE. ) or !
992! chopping ( RND = .FALSE. ) occurs in addition. This may !
993! not be a reliable guide to the way in which the machine !
994! performs its arithmetic. !
995! !
996! EPS (output) DOUBLE PRECISION !
997! The smallest positive number such that !
998! !
999! fl( 1.0 - EPS ) .LT. 1.0, !
1000! !
1001! where fl denotes the computed value. !
1002! !
1003! EMIN (output) INTEGER !
1004! The minimum exponent before (gradual) underflow occurs. !
1005! !
1006! RMIN (output) DOUBLE PRECISION !
1007! The smallest normalized number for the machine, given by !
1008! BASE**( EMIN - 1 ), where BASE is the floating point !
1009! value of BETA. !
1010! !
1011! EMAX (output) INTEGER !
1012! The maximum exponent before overflow occurs. !
1013! !
1014! RMAX (output) DOUBLE PRECISION !
1015! The largest positive number for the machine, given by !
1016! BASE**EMAX * ( 1 - EPS ), where BASE is the floating !
1017! point value of BETA. !
1018! !
1019! Further Details: !
1020! !
1021! The computation of EPS is based on a routine PARANOIA by !
1022! W. Kahan of the University of California at Berkeley. !
1023! !
1024! -- LAPACK auxiliary routine (version 2.0) -- !
1025! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., !
1026! Courant Institute, Argonne National Lab, and Rice University !
1027! October 31, 1992 !
1028! !
1029!=======================================================================
1030!
1031! Imported variable declarations.
1032!
1033 logical, intent(out) :: rnd
1034!
1035 integer, intent(out) :: beta, emax, emin, t
1036!
1037 real (dp), intent(out) :: eps, rmax, rmin
1038!
1039! Local variable declarations.
1040!
1041 logical, save :: first = .true.
1042 logical, save :: iwarn = .false.
1043 logical :: ieee, lieee1, lrnd
1044!
1045 integer, save :: lbeta, lemax, lemin, lt
1046 integer :: gnmin, gpmin, i, ngnmin, ngpmin
1047!
1048 real (dp), save :: leps, lrmax, lrmin
1049 real (dp) :: a, b, c, half, one, rbase, &
1050 & SIXTH, SMALL, THIRD, TWO, ZERO
1051!
1052!-----------------------------------------------------------------------
1053! Executable Statements.
1054!-----------------------------------------------------------------------
1055!
1056 first_pass : IF (first) THEN
1057 first = .false.
1058 zero = 0
1059 one = 1
1060 two = 2
1061!
1062! LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of
1063! BETA, T, RND, EPS, EMIN and RMIN.
1064!
1065! Throughout this routine we use the function DLAMC3 to ensure
1066! that relevant values are stored and not held in registers, or
1067! are not affected by optimizers.
1068!
1069! DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1.
1070!
1071 CALL dlamc1 (lbeta, lt, lrnd, lieee1)
1072!
1073! Start to find EPS.
1074!
1075 b = lbeta
1076 a = b**(-lt)
1077 leps = a
1078!
1079! Try some tricks to see whether or not this is the correct EPS.
1080!
1081 b = two / 3
1082 half = one / 2
1083 sixth = dlamc3(b, -half)
1084 third = dlamc3(sixth, sixth)
1085 b = dlamc3(third, -half)
1086 b = dlamc3(b, sixth)
1087 b = abs(b)
1088 IF (b .LT. leps) b = leps
1089!
1090 leps = 1
1091
1092 DO WHILE ((leps.GT.b) .AND. (b.GT.zero))
1093 leps = b
1094 c = dlamc3(half*leps, (two**5)*(leps**2))
1095 c = dlamc3(half, -c)
1096 b = dlamc3(half, c)
1097 c = dlamc3(half, -b)
1098 b = dlamc3(half, c)
1099 END DO
1100!
1101 IF (a .LT. leps) leps = a
1102!
1103! Computation of EPS complete.
1104!
1105! Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)).
1106! Keep dividing A by BETA until (gradual) underflow occurs. This
1107! is detected when we cannot recover the previous A.
1108!
1109 rbase = one / lbeta
1110 small = one
1111 DO i = 1, 3
1112 small = dlamc3(small*rbase, zero)
1113 END DO
1114 a = dlamc3(one, small)
1115 CALL dlamc4(ngpmin, one, lbeta)
1116 CALL dlamc4(ngnmin, -one, lbeta)
1117 CALL dlamc4(gpmin, a, lbeta)
1118 CALL dlamc4(gnmin, -a, lbeta)
1119 ieee = .false.
1120!
1121 IF ((ngpmin.EQ.ngnmin) .AND. (gpmin.EQ.gnmin)) THEN
1122 IF (ngpmin .EQ. gpmin) THEN
1123 lemin = ngpmin
1124! ( Non twos-complement machines, no gradual underflow;
1125! e.g., VAX )
1126 ELSE IF ((gpmin-ngpmin) .EQ. 3) THEN
1127 lemin = ngpmin - 1 + lt
1128 ieee = .true.
1129! ( Non twos-complement machines, with gradual underflow;
1130! e.g., IEEE standard followers )
1131 ELSE
1132 lemin = min(ngpmin, gpmin)
1133! ( A guess; no known machine )
1134 iwarn = .true.
1135 END IF
1136!
1137 ELSE IF ((ngpmin.EQ.gpmin) .AND. (ngnmin.EQ.gnmin)) THEN
1138 IF (abs(ngpmin-ngnmin) .EQ. 1) THEN
1139 lemin = max(ngpmin, ngnmin)
1140! ( Twos-complement machines, no gradual underflow;
1141! e.g., CYBER 205 )
1142 ELSE
1143 lemin = min(ngpmin, ngnmin)
1144! ( A guess; no known machine )
1145 iwarn = .true.
1146 END IF
1147!
1148 ELSE IF ((abs(ngpmin-ngnmin).EQ.1) .AND. (gpmin.EQ.gnmin)) THEN
1149 IF ((gpmin-min(ngpmin, ngnmin)) .EQ. 3) THEN
1150 lemin = max(ngpmin, ngnmin ) - 1 + lt
1151! ( Twos-complement machines with gradual underflow;
1152! no known machine )
1153 ELSE
1154 lemin = min(ngpmin, ngnmin)
1155! ( A guess; no known machine )
1156 iwarn = .true.
1157 END IF
1158!
1159 ELSE
1160 lemin = min(ngpmin, ngnmin, gpmin, gnmin)
1161! ( A guess; no known machine )
1162 iwarn = .true.
1163 END IF
1164!
1165! Comment out this if block if EMIN is okay.
1166!
1167 IF (iwarn) THEN
1168 first = .true.
1169 print 10, lemin
1170 END IF
1171!
1172! Assume IEEE arithmetic if we found denormalised numbers above,
1173! or if arithmetic seems to round in the IEEE style, determined
1174! in routine DLAMC1. A true IEEE machine should have both things
1175! true; however, faulty machines may have one or the other.
1176!
1177 ieee = ieee .OR. lieee1
1178!
1179! Compute RMIN by successive division by BETA. We could compute
1180! RMIN as BASE**( EMIN - 1 ), but some machines underflow during
1181! this computation.
1182!
1183 lrmin = 1
1184 DO i = 1, 1 - lemin
1185 lrmin = dlamc3(lrmin*rbase, zero)
1186 END DO
1187!
1188! Finally, call DLAMC5 to compute EMAX and RMAX.
1189!
1190 CALL dlamc5 (lbeta, lt, lemin, ieee, lemax, lrmax)
1191 END IF first_pass
1192!
1193 beta = lbeta
1194 t = lt
1195 rnd = lrnd
1196 eps = leps
1197 emin = lemin
1198 rmin = lrmin
1199 emax = lemax
1200 rmax = lrmax
1201!
1202 10 FORMAT ( / / ' WARNING. The value EMIN may be incorrect:-', &
1203 & ' EMIN = ', i0, /, &
1204 & ' If, after inspection, the value EMIN looks', &
1205 & ' acceptable please comment out ', &
1206 & /, ' the IF block as marked within the code of routine', &
1207 & ' DLAMC2,', / ' otherwise supply EMIN explicitly.', / )
1208!
1209 RETURN
1210 END SUBROUTINE dlamc2
1211!
1212 FUNCTION dlamc3 (A, B) RESULT (AplusB)
1213!
1214!=======================================================================
1215! !
1216! DLAMC3 is intended to force A and B to be stored prior to doing !
1217! the addition of A and B , for use in situations where optimizers !
1218! might hold one of these in a register. !
1219! !
1220! Arguments: !
1221! !
1222! A, B (input) DOUBLE PRECISION !
1223! The values A and B. !
1224! !
1225! -- LAPACK auxiliary routine (version 2.0) -- !
1226! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., !
1227! Courant Institute, Argonne National Lab, and Rice University !
1228! October 31, 1992 !
1229! !
1230!=======================================================================
1231!
1232! Imported variable declarations.
1233!
1234 real (dp), intent(in) :: a, b
1235!
1236! Local variable declarations.
1237!
1238 real (dp) :: aplusb
1239!
1240!-----------------------------------------------------------------------
1241! Executable Statements.
1242!-----------------------------------------------------------------------
1243!
1244 aplusb = a + b
1245!
1246 RETURN
1247 END FUNCTION dlamc3
1248!
1249 SUBROUTINE dlamc4 (EMIN, START, BASE)
1250!
1251!=======================================================================
1252! !
1253! DLAMC4 is a service routine for DLAMC2. !
1254! !
1255! Arguments: !
1256! !
1257! EMIN (output) EMIN !
1258! The minimum exponent before (gradual) underflow, computed !
1259! by setting A = START and dividing by BASE until the !
1260! previous A can not be recovered. !
1261! !
1262! START (input) DOUBLE PRECISION !
1263! The starting point for determining EMIN. !
1264! !
1265! BASE (input) INTEGER !
1266! The base of the machine. !
1267! !
1268! -- LAPACK auxiliary routine (version 2.0) -- !
1269! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., !
1270! Courant Institute, Argonne National Lab, and Rice University !
1271! October 31, 1992 !
1272! !
1273!=======================================================================
1274!
1275! Imported variable declarations.
1276!
1277 integer, intent(in ) :: base
1278 integer, intent(out) :: emin
1279!
1280 real (dp), intent(in) :: start
1281!
1282! Local variable declarations.
1283!
1284 integer :: i
1285!
1286 real (dp) :: a, b1, b2, c1, c2, d1, d2, one, rbase, zero
1287!
1288!-----------------------------------------------------------------------
1289! Executable Statements.
1290!-----------------------------------------------------------------------
1291!
1292 a = start
1293 one = 1
1294 rbase = one / base
1295 zero = 0
1296 emin = 1
1297 b1 = dlamc3(a*rbase, zero)
1298 c1 = a
1299 c2 = a
1300 d1 = a
1301 d2 = a
1302!
1303 DO WHILE ((c1.EQ.a) .AND. (c2.EQ.a) .AND. &
1304 & (d1.EQ.a) .AND. (d2.EQ.a))
1305 emin = emin - 1
1306 a = b1
1307 b1 = dlamc3(a / base, zero)
1308 c1 = dlamc3(b1*base, zero)
1309 d1 = zero
1310 DO i = 1, base
1311 d1 = d1 + b1
1312 END DO
1313 b2 = dlamc3(a*rbase, zero)
1314 c2 = dlamc3(b2 / rbase, zero)
1315 d2 = zero
1316 DO i = 1, base
1317 d2 = d2 + b2
1318 END DO
1319 END DO
1320!
1321 RETURN
1322 END SUBROUTINE dlamc4
1323!
1324 SUBROUTINE dlamc5 (BETA, P, EMIN, IEEE, EMAX, RMAX)
1325!
1326!=======================================================================
1327! !
1328! DLAMC5 attempts to compute RMAX, the largest machine floating-point !
1329! number, without overflow. It assumes that EMAX + abs(EMIN) sum !
1330! approximately to a power of 2. It will fail on machines where this !
1331! assumption does not hold, for example, the Cyber 205 (EMIN = -28625,!
1332! EMAX = 28718). It will also fail if the value supplied for EMIN is !
1333! too large (i.e. too close to zero), probably with overflow. !
1334! !
1335! Arguments: !
1336! !
1337! BETA (input) INTEGER !
1338! The base of floating-point arithmetic. !
1339! !
1340! P (input) INTEGER !
1341! The number of base BETA digits in the mantissa of a !
1342! floating-point value. !
1343! !
1344! EMIN (input) INTEGER !
1345! The minimum exponent before (gradual) underflow. !
1346! !
1347! IEEE (input) LOGICAL !
1348! A logical flag specifying whether or not the arithmetic !
1349! system is thought to comply with the IEEE standard. !
1350! !
1351! EMAX (output) INTEGER !
1352! The largest exponent before overflow !
1353! !
1354! RMAX (output) DOUBLE PRECISION !
1355! The largest machine floating-point number. !
1356! !
1357! -- LAPACK auxiliary routine (version 2.0) -- !
1358! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., !
1359! Courant Institute, Argonne National Lab, and Rice University !
1360! October 31, 1992 !
1361! !
1362!=======================================================================
1363!
1364! Imported variable declarations.
1365!
1366 logical, intent(in) :: ieee
1367!
1368 integer, intent(in) :: beta, emin, p
1369 integer, intent(out) :: emax
1370!
1371 real (dp), intent(out) :: rmax
1372!
1373! Local variable declarations.
1374!
1375 real (dp), parameter :: zero = 0.0_dp
1376 real (dp), parameter :: one = 1.0_dp
1377!
1378 integer :: exbits, expsum, i, lexp, nbits, try, uexp
1379!
1380 real (dp) :: oldy, recbas, y, z
1381!
1382!-----------------------------------------------------------------------
1383! Executable Statements.
1384!-----------------------------------------------------------------------
1385!
1386! First compute LEXP and UEXP, two powers of 2 that bound
1387! ABS(EMIN). We then assume that EMAX + ABS(EMIN) will sum
1388! approximately to the bound that is closest to ABS(EMIN).
1389! (EMAX is the exponent of the required number RMAX).
1390!
1391 lexp = 1
1392 exbits = 1
1393#ifdef REMOVE_LAPACK_GOTOS
1394 try_loop : DO ! forever DO-loop
1395 try = lexp*2
1396 IF (try .LE. (-emin)) THEN
1397 lexp = try
1398 exbits = exbits + 1
1399 cycle try_loop ! iterate DO-loop again
1400 ELSE
1401 EXIT try_loop ! terminate DO-loop
1402 END IF
1403 END DO try_loop
1404#else
1405 10 CONTINUE
1406 try = lexp*2
1407 IF (try .LE. (-emin)) THEN
1408 lexp = try
1409 exbits = exbits + 1
1410 GO TO 10
1411 END IF
1412#endif
1413 IF (lexp .EQ. -emin) THEN
1414 uexp = lexp
1415 ELSE
1416 uexp = try
1417 exbits = exbits + 1
1418 END IF
1419!
1420! Now -LEXP is less than or equal to EMIN, and -UEXP is greater
1421! than or equal to EMIN. EXBITS is the number of bits needed to
1422! store the exponent.
1423!
1424 IF ((uexp+emin) .GT. (-lexp-emin)) THEN
1425 expsum = 2*lexp
1426 ELSE
1427 expsum = 2*uexp
1428 END IF
1429!
1430! EXPSUM is the exponent range, approximately equal to
1431! EMAX - EMIN + 1 .
1432!
1433 emax = expsum + emin - 1
1434 nbits = 1 + exbits + p
1435!
1436! NBITS is the total number of bits needed to store a
1437! floating-point number.
1438!
1439 IF ((mod(nbits, 2).EQ.1) .AND. (beta.EQ.2)) THEN
1440!
1441! Either there are an odd number of bits used to store a
1442! floating-point number, which is unlikely, or some bits are
1443! not used in the representation of numbers, which is possible,
1444! (e.g. Cray machines) or the mantissa has an implicit bit,
1445! (e.g. IEEE machines, Dec Vax machines), which is perhaps the
1446! most likely. We have to assume the last alternative.
1447! If this is true, then we need to reduce EMAX by one because
1448! there must be some way of representing zero in an implicit-bit
1449! system. On machines like Cray, we are reducing EMAX by one
1450! unnecessarily.
1451!
1452 emax = emax - 1
1453 END IF
1454!
1455 IF (ieee) THEN
1456!
1457! Assume we are on an IEEE machine which reserves one exponent
1458! for infinity and NaN.
1459!
1460 emax = emax - 1
1461 END IF
1462!
1463! Now create RMAX, the largest machine number, which should
1464! be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
1465!
1466! First compute 1.0 - BETA**(-P), being careful that the
1467! result is less than 1.0 .
1468!
1469 recbas = one / beta
1470 z = beta - one
1471 y = zero
1472 DO i = 1, p
1473 z = z*recbas
1474 IF (y .LT. one) oldy = y
1475 y = dlamc3(y, z)
1476 END DO
1477 IF (y .GE. one) y = oldy
1478!
1479! Now multiply by BETA**EMAX to get RMAX.
1480!
1481 DO i = 1, emax
1482 y = dlamc3(y*beta, zero)
1483 END DO
1484!
1485 rmax = y
1486!
1487 RETURN
1488 END SUBROUTINE dlamc5
1489!
1490 FUNCTION dlamch (CMACH) RESULT (RMACH)
1491!
1492!=======================================================================
1493! !
1494! DLAMCH determines double precision machine parameters. !
1495! !
1496! Arguments: !
1497! !
1498! CMACH (input) CHARACTER*1 !
1499! Specifies the value to be returned by DLAMCH: !
1500! = 'E' or 'e', DLAMCH := eps !
1501! = 'S' or 's , DLAMCH := sfmin !
1502! = 'B' or 'b', DLAMCH := base !
1503! = 'P' or 'p', DLAMCH := eps*base !
1504! = 'N' or 'n', DLAMCH := t !
1505! = 'R' or 'r', DLAMCH := rnd !
1506! = 'M' or 'm', DLAMCH := emin !
1507! = 'U' or 'u', DLAMCH := rmin !
1508! = 'L' or 'l', DLAMCH := emax !
1509! = 'O' or 'o', DLAMCH := rmax !
1510! !
1511! where !
1512! !
1513! eps = relative machine precision !
1514! sfmin = safe minimum, such that 1/sfmin does not overflow !
1515! base = base of the machine !
1516! prec = eps*base !
1517! t = number of (base) digits in the mantissa !
1518! rnd = 1.0 when rounding occurs in addition, 0.0 otherwise !
1519! emin = minimum exponent before (gradual) underflow !
1520! rmin = underflow threshold - base**(emin-1) !
1521! emax = largest exponent before overflow !
1522! rmax = overflow threshold - (base**emax)*(1-eps) !
1523! !
1524! -- LAPACK auxiliary routine (version 2.0) -- !
1525! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., !
1526! Courant Institute, Argonne National Lab, and Rice University !
1527! October 31, 1992 !
1528! !
1529!=======================================================================
1530!
1531! Imported variable declarations.
1532!
1533 character (len=1), intent(in) :: cmach
1534!
1535! Local variable declarations.
1536!
1537 real (dp), parameter :: zero = 0.0_dp
1538 real (dp), parameter :: one = 1.0_dp
1539!
1540 logical, save :: first = .true.
1541 logical :: lrnd
1542!
1543 integer :: beta, imax, imin, it
1544!
1545 real (dp), save :: eps, sfmin, base, t, rnd, emin, rmin, &
1546 & EMAX, RMAX, PREC
1547 real (dp) :: rmach, small
1548!
1549!-----------------------------------------------------------------------
1550! Executable Statements ..
1551!-----------------------------------------------------------------------
1552!
1553 IF (first) THEN
1554 first = .false.
1555 CALL dlamc2 (beta, it, lrnd, eps, imin, rmin, imax, rmax)
1556 base = beta
1557 t = it
1558 IF (lrnd) THEN
1559 rnd = one
1560 eps = (base**(1-it)) / 2
1561 ELSE
1562 rnd = zero
1563 eps = base**(1-it)
1564 END IF
1565 prec = eps*base
1566 emin = imin
1567 emax = imax
1568 sfmin = rmin
1569 small = one / rmax
1570 IF (small.GE.sfmin) THEN
1571!
1572! Use SMALL plus a bit, to avoid the possibility of rounding
1573! causing overflow when computing 1/sfmin.
1574!
1575 sfmin = small*(one+eps)
1576 END IF
1577 END IF
1578!
1579 IF (lsame(cmach, 'E')) THEN
1580 rmach = eps
1581 ELSE IF (lsame(cmach, 'S')) THEN
1582 rmach = sfmin
1583 ELSE IF (lsame(cmach, 'B')) THEN
1584 rmach = base
1585 ELSE IF (lsame(cmach, 'P')) THEN
1586 rmach = prec
1587 ELSE IF (lsame(cmach, 'N')) THEN
1588 rmach = t
1589 ELSE IF (lsame(cmach, 'R')) THEN
1590 rmach = rnd
1591 ELSE IF (lsame(cmach, 'M')) THEN
1592 rmach = emin
1593 ELSE IF (lsame(cmach, 'U')) THEN
1594 rmach = rmin
1595 ELSE IF (lsame(cmach, 'L')) THEN
1596 rmach = emax
1597 ELSE IF (lsame(cmach, 'O')) THEN
1598 rmach = rmax
1599 END IF
1600!
1601 RETURN
1602 END FUNCTION dlamch
1603!
1604 FUNCTION dlanst (NORM, N, D, E)
1605!
1606!=======================================================================
1607! !
1608! DLANST returns the value of the one norm, or the Frobenius norm, !
1609! or the infinity norm, or the element of largest absolute value of !
1610! a real symmetric tridiagonal matrix A. !
1611! !
1612! Description: !
1613! !
1614! DLANST returns the value !
1615! !
1616! ANORM = ( MAX(ABS(A(i,j))), NORM = 'M' or 'm' !
1617! ( !
1618! ( norm1(A), NORM = '1', 'O' or 'o' !
1619! ( !
1620! ( normI(A), NORM = 'I' or 'i' !
1621! ( !
1622! ( normF(A), NORM = 'F', 'f', 'E' or 'e' !
1623! where !
1624! !
1625! norm1 denotes the one norm of a matrix (maximum column sum), !
1626! normI denotes the infinity norm of a matrix (maximum row sum), and !
1627! normF denotes the Frobenius norm of a matrix (square root of sum of !
1628! squares). !
1629! Note that MAX(ABS(A(i,j))) is not a matrix norm. !
1630! !
1631! Arguments: !
1632! !
1633! NORM (input) CHARACTER*1 !
1634! Specifies the value to be returned in DLANST as described !
1635! above. !
1636! !
1637! N (input) INTEGER !
1638! The order of the matrix A. N >= 0. When N = 0, DLANST is !
1639! set to zero. !
1640! !
1641! D (input) DOUBLE PRECISION array, dimension (N) !
1642! The diagonal elements of A. !
1643! !
1644! E (input) DOUBLE PRECISION array, dimension (N-1) !
1645! The (n-1) sub-diagonal or super-diagonal elements of A. !
1646! !
1647! -- LAPACK auxiliary routine (version 2.0) -- !
1648! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., !
1649! Courant Institute, Argonne National Lab, and Rice University !
1650! February 29, 1992 !
1651! !
1652!=======================================================================
1653!
1654! Imported variable declarations.
1655!
1656 integer, intent(in) :: n
1657!
1658 real (dp), intent(in) :: d(*), e(*)
1659!
1660 character (len=1), intent(in) :: norm
1661!
1662! Local variable declarations.
1663!
1664 real (dp), parameter :: zero = 0.0_dp
1665 real (dp), parameter :: one = 1.0_dp
1666!
1667 integer :: i
1668 real (dp) :: anorm, scale, sum
1669!
1670!-----------------------------------------------------------------------
1671! Executable Statements.
1672!-----------------------------------------------------------------------
1673!
1674 IF (n.LE.0) THEN
1675 anorm = zero
1676 ELSE IF (lsame(norm, 'M')) THEN
1677!
1678! Find MAX(ABS(A(i,j))).
1679!
1680 anorm = abs(d(n))
1681 DO i = 1, n - 1
1682 anorm = max(anorm, abs(d(i)))
1683 anorm = max(anorm, abs(e(i)))
1684 END DO
1685 ELSE IF (lsame(norm, 'O') .OR. norm.EQ.'1' .OR. &
1686 & lsame(norm, 'I')) THEN
1687!
1688! Find norm1(A).
1689!
1690 IF (n.EQ.1) THEN
1691 anorm = abs(d(1))
1692 ELSE
1693 anorm = max(abs(d(1)) + abs(e(1)), &
1694 & abs(e(n-1)) + abs(d(n)))
1695 DO i = 2, n - 1
1696 anorm = max(anorm, abs(d(i)) + abs(e(i)) + abs(e(i-1)))
1697 END DO
1698 END IF
1699 ELSE IF ((lsame(norm, 'F')) .OR. (lsame(norm, 'E'))) THEN
1700!
1701! Find normF(A).
1702!
1703 scale = zero
1704 sum = one
1705 IF (n .GT. 1) THEN
1706 CALL dlassq (n-1, e, 1, scale, sum)
1707 sum = 2*sum
1708 END IF
1709 CALL dlassq (n, d, 1, scale, sum)
1710 anorm = scale*sqrt(sum)
1711 END IF
1712!
1713 dlanst = anorm
1714!
1715 RETURN
1716 END FUNCTION dlanst
1717!
1718 FUNCTION dlapy2 (X, Y) RESULT (LAPY2)
1719!
1720!=======================================================================
1721! !
1722! DLAPY2 returns SQRT(x**2+y**2), taking care not to cause !
1723! unnecessary overflow. !
1724! !
1725! Arguments: !
1726! !
1727! X (input) DOUBLE PRECISION !
1728! Y (input) DOUBLE PRECISION !
1729! X and Y specify the values x and y. !
1730! !
1731! -- LAPACK auxiliary routine (version 2.0) -- !
1732! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., !
1733! Courant Institute, Argonne National Lab, and Rice University !
1734! October 31, 1992 !
1735! !
1736!=======================================================================
1737!
1738! Imported variable declarations.
1739!
1740 real (dp), intent(in) :: x, y
1741!
1742! Local variable declarations.
1743!
1744 real (dp), parameter :: zero = 0.0_dp
1745 real (dp), parameter :: one = 1.0_dp
1746!
1747 real (dp) :: w, xabs, yabs, z
1748 real (dp) :: lapy2
1749!
1750!-----------------------------------------------------------------------
1751! Executable Statements ..
1752!-----------------------------------------------------------------------
1753!
1754 xabs = abs(x)
1755 yabs = abs(y)
1756 w = max(xabs, yabs)
1757 z = min(xabs, yabs)
1758 IF (z .EQ. zero) THEN
1759 lapy2 = w
1760 ELSE
1761 lapy2 = w*sqrt(one+(z / w)**2)
1762 END IF
1763!
1764 RETURN
1765 END FUNCTION dlapy2
1766!
1767 SUBROUTINE dlartg (F, G, CS, SN, R)
1768!
1769!=======================================================================
1770! !
1771! DLARTG generate a plane rotation so that !
1772! !
1773! [ CS SN ] . [ F ] = [ R ] where CS**2 + SN**2 = 1. !
1774! [ -SN CS ] [ G ] [ 0 ] !
1775! !
1776! This is a slower, more accurate version of the BLAS1 routine DROTG, !
1777! with the following other differences: !
1778! !
1779! F and G are unchanged on return. !
1780! If G=0, then CS=1 and SN=0. !
1781! If F=0 and (G .ne. 0), then CS=0 and SN=1 without doing any !
1782! floating point operations (saves work in DBDSQR when !
1783! there are zeros on the diagonal). !
1784! !
1785! If F exceeds G in magnitude, CS will be positive. !
1786! !
1787! Arguments: !
1788! !
1789! F (input) DOUBLE PRECISION !
1790! The first component of vector to be rotated. !
1791! !
1792! G (input) DOUBLE PRECISION !
1793! The second component of vector to be rotated. !
1794! !
1795! CS (output) DOUBLE PRECISION !
1796! The cosine of the rotation. !
1797! !
1798! SN (output) DOUBLE PRECISION !
1799! The sine of the rotation. !
1800! !
1801! R (output) DOUBLE PRECISION !
1802! The nonzero component of the rotated vector. !
1803! !
1804! -- LAPACK auxiliary routine (version 2.0) -- !
1805! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., !
1806! Courant Institute, Argonne National Lab, and Rice University !
1807! September 30, 1994 !
1808! !
1809!=======================================================================
1810!
1811! Imported variable declarations.
1812!
1813 real (dp), intent(in ) :: f, g
1814 real (dp), intent(out) :: cs, r, sn
1815!
1816! Local variable declarations.
1817!
1818 real (dp), parameter :: zero = 0.0_dp
1819 real (dp), parameter :: one = 1.0_dp
1820 real (dp), parameter :: two = 2.0_dp
1821!
1822 logical, save :: first = .true.
1823!
1824 integer :: count, i
1825!
1826 real (dp), save :: safmx2, safmin, safmn2
1827 real (dp) :: eps, f1, g1, scale
1828!
1829!-----------------------------------------------------------------------
1830! Executable Statements.
1831!-----------------------------------------------------------------------
1832!
1833 IF (first) THEN
1834 first = .false.
1835 safmin = dlamch('S')
1836 eps = dlamch('E')
1837 safmn2 = dlamch('B')**int(log(safmin / eps) / &
1838 & log(dlamch('B')) / two)
1839 safmx2 = one / safmn2
1840 END IF
1841 IF (g .EQ. zero) THEN
1842 cs = one
1843 sn = zero
1844 r = f
1845 ELSE IF (f .EQ. zero) THEN
1846 cs = zero
1847 sn = one
1848 r = g
1849 ELSE
1850 f1 = f
1851 g1 = g
1852 scale = max(abs(f1), abs(g1))
1853 IF (scale .GE. safmx2) THEN
1854 count = 0
1855#ifdef REMOVE_LAPACK_GOTOS
1856 loop_1 : DO ! forever DO-loop
1857 count = count + 1
1858 f1 = f1*safmn2
1859 g1 = g1*safmn2
1860 scale = max(abs(f1), abs(g1))
1861 IF (scale .GE. safmx2) cycle loop_1 ! iterate DO-loop again
1862 r = sqrt(f1**2 + g1**2)
1863 cs = f1 / r
1864 sn = g1 / r
1865 DO i = 1, count
1866 r = r*safmx2
1867 END DO
1868 EXIT loop_1 ! terminate DO-loop
1869 END DO loop_1
1870#else
1871 10 CONTINUE
1872 count = count + 1
1873 f1 = f1*safmn2
1874 g1 = g1*safmn2
1875 scale = max(abs(f1), abs(g1))
1876 IF (scale .GE. safmx2) GO TO 10
1877 r = sqrt(f1**2 + g1**2)
1878 cs = f1 / r
1879 sn = g1 / r
1880 DO i = 1, count
1881 r = r*safmx2
1882 END DO
1883#endif
1884 ELSE IF (scale .LE. safmn2) THEN
1885 count = 0
1886#ifdef REMOVE_LAPACK_GOTOS
1887 loop_2 : DO ! forever DO-loop
1888 count = count + 1
1889 f1 = f1*safmx2
1890 g1 = g1*safmx2
1891 scale = max(abs(f1), abs(g1))
1892 IF (scale .LE. safmn2) cycle loop_2 ! iterate DO-loop again
1893 r = sqrt(f1**2 + g1**2)
1894 cs = f1 / r
1895 sn = g1 / r
1896 DO i = 1, count
1897 r = r*safmn2
1898 END DO
1899 EXIT loop_2 ! terminate DO-loop
1900 END DO loop_2
1901#else
1902 30 CONTINUE
1903 count = count + 1
1904 f1 = f1*safmx2
1905 g1 = g1*safmx2
1906 scale = max(abs(f1), abs(g1))
1907 IF (scale .LE. safmn2) GO TO 30
1908 r = sqrt(f1**2 + g1**2)
1909 cs = f1 / r
1910 sn = g1 / r
1911 DO i = 1, count
1912 r = r*safmn2
1913 END DO
1914#endif
1915 ELSE
1916 r = sqrt(f1**2 + g1**2)
1917 cs = f1 / r
1918 sn = g1 / r
1919 END IF
1920 IF ((abs(f).GT.abs(g)) .AND. (cs.LT.zero)) THEN
1921 cs = -cs
1922 sn = -sn
1923 r = -r
1924 END IF
1925 END IF
1926!
1927 RETURN
1928 END SUBROUTINE dlartg
1929!
1930 SUBROUTINE dlascl (TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
1931!
1932!=======================================================================
1933! !
1934! DLASCL multiplies the M by N real matrix A by the real scalar !
1935! CTO/CFROM. This is done without over/underflow as long as the !
1936! final result CTO*A(I,J)/CFROM does not over/underflow. TYPE !
1937! specifies that A may be full, upper triangular, lower triangular, !
1938! upper Hessenberg, or banded. !
1939! !
1940! Arguments: !
1941! !
1942! TYPE (input) CHARACTER*1 !
1943! TYPE indices the storage type of the input matrix. !
1944! = 'G': A is a full matrix. !
1945! = 'L': A is a lower triangular matrix. !
1946! = 'U': A is an upper triangular matrix. !
1947! = 'H': A is an upper Hessenberg matrix. !
1948! = 'B': A is a symmetric band matrix with lower bandwidth !
1949! KL and upper bandwidth KU and with the only the !
1950! lower half stored. !
1951! = 'Q': A is a symmetric band matrix with lower bandwidth !
1952! KL and upper bandwidth KU and with the only the !
1953! upper half stored. !
1954! = 'Z': A is a band matrix with lower bandwidth KL and !
1955! upper bandwidth KU. !
1956! !
1957! KL (input) INTEGER !
1958! The lower bandwidth of A. Referenced only if TYPE = 'B', !
1959! 'Q' or 'Z'. !
1960! !
1961! KU (input) INTEGER !
1962! The upper bandwidth of A. Referenced only if TYPE = 'B', !
1963! 'Q' or 'Z'. !
1964! !
1965! CFROM (input) DOUBLE PRECISION !
1966! CTO (input) DOUBLE PRECISION !
1967! The matrix A is multiplied by CTO/CFROM. A(I,J) is computed !
1968! without over/underflow if the final result CTO*A(I,J)/CFROM !
1969! can be represented without over/underflow. CFROM must be !
1970! nonzero. !
1971! !
1972! M (input) INTEGER !
1973! The number of rows of the matrix A. M >= 0. !
1974! !
1975! N (input) INTEGER !
1976! The number of columns of the matrix A. N >= 0. !
1977! !
1978! A (input/output) DOUBLE PRECISION array, dimension (LDA,M) !
1979! The matrix to be multiplied by CTO/CFROM. See TYPE for the !
1980! storage type. !
1981! !
1982! LDA (input) INTEGER !
1983! The leading dimension of the array A. LDA >= max(1,M). !
1984! !
1985! INFO (output) INTEGER !
1986! 0 - successful exit !
1987! <0 - if INFO = -i, the i-th argument had an illegal value. !
1988! !
1989! -- LAPACK auxiliary routine (version 2.0) -- !
1990! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., !
1991! Courant Institute, Argonne National Lab, and Rice University !
1992! February 29, 1992 !
1993! !
1994!=======================================================================
1995!
1996! Imported variable declarations.
1997!
1998!
1999 integer, intent (in ) :: kl, ku, lda, m, n
2000 integer, intent (out) :: info
2001!
2002 real (dp), intent(in ) :: cfrom, cto
2003 real (dp), intent(inout) :: a(lda, *)
2004!
2005 character (len=1), intent(in) :: type
2006!
2007! Local variable declarations.
2008!
2009 real (r8), parameter :: zero = 0.0_dp
2010 real (r8), parameter :: one = 1.0_dp
2011!
2012 logical :: done
2013!
2014 integer :: i, itype, j, k1, k2, k3, k4
2015!
2016 real (dp) :: bignum, cfrom1, cfromc, cto1, ctoc, mul, smlnum
2017!
2018!-----------------------------------------------------------------------
2019! Executable Statements.
2020!-----------------------------------------------------------------------
2021!
2022! Test the input arguments.
2023!
2024 info = 0
2025!
2026 IF (lsame(TYPE, 'G')) then
2027 itype = 0
2028 ELSE IF (lsame(TYPE, 'L')) then
2029 itype = 1
2030 ELSE IF (lsame(TYPE, 'U')) then
2031 itype = 2
2032 ELSE IF (lsame(TYPE, 'H')) then
2033 itype = 3
2034 ELSE IF (lsame(TYPE, 'B')) then
2035 itype = 4
2036 ELSE IF (lsame(TYPE, 'Q')) then
2037 itype = 5
2038 ELSE IF (lsame(TYPE, 'Z')) then
2039 itype = 6
2040 ELSE
2041 itype = -1
2042 END IF
2043!
2044 IF (itype .EQ. -1) THEN
2045 info = -1
2046 ELSE IF (cfrom .EQ. zero) THEN
2047 info = -4
2048 ELSE IF (m .LT. 0) THEN
2049 info = -6
2050 ELSE IF (n.LT.0 .OR. (itype.EQ.4 .AND. n.NE.m) .OR. &
2051 & (itype.EQ.5 .AND. n.NE.m)) THEN
2052 info = -7
2053 ELSE IF (itype.LE.3 .AND. lda.LT.max(1, m)) THEN
2054 info = -9
2055 ELSE IF (itype .GE. 4) THEN
2056 IF (kl.LT.0 .OR. kl.GT.max(m-1, 0)) THEN
2057 info = -2
2058 ELSE IF (ku.LT.0 .OR. ku.GT.max(n-1, 0) .OR. &
2059 & ((itype.EQ.4 .OR. itype.EQ.5) .AND. kl.NE.ku)) THEN
2060 info = -3
2061 ELSE IF ((itype.EQ.4 .AND. lda.LT.kl+1) .OR. &
2062 & (itype.EQ.5 .AND. lda.LT.ku+1 ) .OR. &
2063 & (itype.EQ.6 .AND. lda.LT.2*kl+ku+1)) THEN
2064 info = -9
2065 END IF
2066 END IF
2067!
2068 IF (info .NE. 0) THEN
2069 CALL xerbla ('DLASCL', -info)
2070 RETURN
2071 END IF
2072!
2073! Quick return if possible.
2074!
2075 IF ((n.EQ.0) .OR. (m.EQ.0)) RETURN
2076!
2077! Get machine parameters
2078!
2079 smlnum = dlamch('S')
2080 bignum = one / smlnum
2081!
2082 cfromc = cfrom
2083 ctoc = cto
2084!
2085 done = .false.
2086!
2087 iterate : DO WHILE (.NOT. done)
2088 cfrom1 = cfromc*smlnum
2089 cto1 = ctoc / bignum
2090 IF (abs(cfrom1).GT.abs(ctoc) .AND. ctoc.NE.zero) THEN
2091 mul = smlnum
2092 done = .false.
2093 cfromc = cfrom1
2094 ELSE IF (abs(cto1) .GT. abs(cfromc)) THEN
2095 mul = bignum
2096 done = .false.
2097 ctoc = cto1
2098 ELSE
2099 mul = ctoc / cfromc
2100 done = .true.
2101 END IF
2102!
2103 IF (itype .EQ. 0) THEN
2104!
2105! Full matrix.
2106!
2107 DO j = 1, n
2108 DO i = 1, m
2109 a(i,j) = a(i,j)*mul
2110 END DO
2111 END DO
2112!
2113 ELSE IF (itype .EQ. 1) THEN
2114!
2115! Lower triangular matrix.
2116!
2117 DO j = 1, n
2118 DO i = j, m
2119 a(i,j) = a(i,j)*mul
2120 END DO
2121 END DO
2122!
2123 ELSE IF (itype.EQ.2) THEN
2124!
2125! Upper triangular matrix.
2126!
2127 DO j = 1, n
2128 DO i = 1, min(j, m)
2129 a(i,j) = a(i,j)*mul
2130 END DO
2131 END DO
2132!
2133 ELSE IF (itype .EQ. 3) THEN
2134!
2135! Upper Hessenberg matrix.
2136!
2137 DO j = 1, n
2138 DO i = 1, min(j+1, m)
2139 a(i,j) = a(i,j)*mul
2140 END DO
2141 END DO
2142!
2143 ELSE IF (itype .EQ. 4) THEN
2144!
2145! Lower half of a symmetric band matrix.
2146!
2147 k3 = kl + 1
2148 k4 = n + 1
2149 DO j = 1, n
2150 DO i = 1, min(k3, k4-j)
2151 a(i,j) = a(i,j)*mul
2152 END DO
2153 END DO
2154!
2155 ELSE IF (itype .EQ. 5) THEN
2156!
2157! Upper half of a symmetric band matrix.
2158!
2159 k1 = ku + 2
2160 k3 = ku + 1
2161 DO j = 1, n
2162 DO i = max(k1-j, 1), k3
2163 a(i,j) = a(i,j)*mul
2164 END DO
2165 END DO
2166!
2167 ELSE IF (itype .EQ. 6) THEN
2168!
2169! Band matrix.
2170!
2171 k1 = kl + ku + 2
2172 k2 = kl + 1
2173 k3 = 2*kl + ku + 1
2174 k4 = kl + ku + 1 + m
2175 DO j = 1, n
2176 DO i = max(k1-j, k2), min(k3, k4-j)
2177 a(i,j) = a(i,j)*mul
2178 END DO
2179 END DO
2180!
2181 END IF
2182 END DO iterate
2183!
2184 RETURN
2185 END SUBROUTINE dlascl
2186!
2187 SUBROUTINE dlaset (UPLO, M, N, ALPHA, BETA, A, LDA)
2188!
2189!=======================================================================
2190! !
2191! DLASET initializes an m-by-n matrix A to BETA on the diagonal and !
2192! ALPHA on the offdiagonals. !
2193! !
2194! Arguments: !
2195! !
2196! UPLO (input) CHARACTER*1 !
2197! Specifies the part of the matrix A to be set. !
2198! = 'U': Upper triangular part is set; the strictly !
2199! lower triangular part of A is not changed. !
2200! = 'L': Lower triangular part is set; the strictly !
2201! upper triangular part of A is not changed. !
2202! Otherwise: All of the matrix A is set. !
2203! !
2204! M (input) INTEGER !
2205! The number of rows of the matrix A. M >= 0. !
2206! !
2207! N (input) INTEGER !
2208! The number of columns of the matrix A. N >= 0. !
2209! !
2210! ALPHA (input) DOUBLE PRECISION !
2211! The constant to which the offdiagonal elements are to be !
2212! set. !
2213! !
2214! BETA (input) DOUBLE PRECISION !
2215! The constant to which the diagonal elements are to be set. !
2216! !
2217! A (input/output) DOUBLE PRECISION array, dimension (LDA,N) !
2218! On exit, the leading m-by-n submatrix of A is set as !
2219! follows: !
2220! !
2221! if UPLO = 'U', A(i,j) = ALPHA, 1<=i<=j-1, 1<=j<=n, !
2222! if UPLO = 'L', A(i,j) = ALPHA, j+1<=i<=m, 1<=j<=n, !
2223! otherwise, A(i,j) = ALPHA, 1<=i<=m, 1<=j<=n, i.ne.j, !
2224! !
2225! and, for all UPLO, A(i,i) = BETA, 1<=i<=min(m,n). !
2226! !
2227! LDA (input) INTEGER !
2228! The leading dimension of the array A. LDA >= max(1,M). !
2229! !
2230! !
2231! -- LAPACK auxiliary routine (version 2.0) -- !
2232! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., !
2233! Courant Institute, Argonne National Lab, and Rice University !
2234! October 31, 1992 !
2235! !
2236!=======================================================================
2237!
2238! Imported variable declarations.
2239!
2240 integer, intent(in) :: lda, m, n
2241!
2242 real (dp), intent(in) :: alpha, beta
2243 real (dp), intent(inout) :: a(lda, *)
2244!
2245 character (len=1), intent(in) :: uplo
2246!
2247! Local variable declarations.
2248!
2249 integer :: i, j
2250!
2251!-----------------------------------------------------------------------
2252! Executable Statements.
2253!-----------------------------------------------------------------------
2254!
2255 IF (lsame(uplo, 'U')) THEN
2256!
2257! Set the strictly upper triangular or trapezoidal part of the
2258! array to ALPHA.
2259!
2260 DO j = 2, n
2261 DO i = 1, min(j-1, m)
2262 a(i,j) = alpha
2263 END DO
2264 END DO
2265!
2266 ELSE IF (lsame(uplo, 'L')) THEN
2267!
2268! Set the strictly lower triangular or trapezoidal part of the
2269! array to ALPHA.
2270!
2271 DO j = 1, min(m, n)
2272 DO i = j + 1, m
2273 a(i,j) = alpha
2274 END DO
2275 END DO
2276!
2277 ELSE
2278!
2279! Set the leading m-by-n submatrix to ALPHA.
2280!
2281 DO j = 1, n
2282 DO i = 1, m
2283 a(i,j) = alpha
2284 END DO
2285 END DO
2286 END IF
2287!
2288! Set the first MIN(M,N) diagonal elements to BETA.
2289!
2290 DO i = 1, min(m, n)
2291 a(i,i) = beta
2292 END DO
2293!
2294 RETURN
2295 END SUBROUTINE dlaset
2296!
2297 SUBROUTINE dlasr (SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
2298!
2299!=======================================================================
2300! !
2301! DLASR performs the transformation: !
2302! !
2303! A := P*A, when SIDE = 'L' or 'l' ( Left-hand side ) !
2304! !
2305! A := A*P', when SIDE = 'R' or 'r' ( Right-hand side ) !
2306! !
2307! where A is an m by n real matrix and P is an orthogonal matrix, !
2308! consisting of a sequence of plane rotations determined by the !
2309! parameters PIVOT and DIRECT as follows ( z = m when SIDE = 'L' !
2310! or 'l' and z = n when SIDE = 'R' or 'r' ): !
2311! !
2312! When DIRECT = 'F' or 'f' ( Forward sequence ) then !
2313! !
2314! P = P( z - 1 )*...*P( 2 )*P( 1 ), !
2315! !
2316! and when DIRECT = 'B' or 'b' ( Backward sequence ) then !
2317! !
2318! P = P( 1 )*P( 2 )*...*P( z - 1 ), !
2319! !
2320! where P( k ) is a plane rotation matrix for the following planes: !
2321! !
2322! when PIVOT = 'V' or 'v' ( Variable pivot ), !
2323! the plane ( k, k + 1 ) !
2324! !
2325! when PIVOT = 'T' or 't' ( Top pivot ), !
2326! the plane ( 1, k + 1 ) !
2327! !
2328! when PIVOT = 'B' or 'b' ( Bottom pivot ), !
2329! the plane ( k, z ) !
2330! !
2331! c( k ) and s( k ) must contain the cosine and sine that define !
2332! the matrix P( k ). The two by two plane rotation part of the !
2333! matrix P( k ), R( k ), is assumed to be of the form !
2334! !
2335! R( k ) = ( c( k ) s( k ) ). !
2336! ( -s( k ) c( k ) ) !
2337! !
2338! This version vectorises across rows of the array A when !
2339! SIDE = 'L'. !
2340! !
2341! Arguments: !
2342! !
2343! SIDE (input) CHARACTER*1 !
2344! Specifies whether the plane rotation matrix P is applied !
2345! to A on the left or the right. !
2346! = 'L': Left, compute A := P*A !
2347! = 'R': Right, compute A:= A*P' !
2348! !
2349! DIRECT (input) CHARACTER*1 !
2350! Specifies whether P is a forward or backward sequence of !
2351! plane rotations. !
2352! = 'F': Forward, P = P( z - 1 )*...*P( 2 )*P( 1 ) !
2353! = 'B': Backward, P = P( 1 )*P( 2 )*...*P( z - 1 ) !
2354! !
2355! PIVOT (input) CHARACTER*1 !
2356! Specifies the plane for which P(k) is a plane rotation !
2357! matrix. !
2358! = 'V': Variable pivot, the plane (k,k+1) !
2359! = 'T': Top pivot, the plane (1,k+1) !
2360! = 'B': Bottom pivot, the plane (k,z) !
2361! !
2362! M (input) INTEGER !
2363! The number of rows of the matrix A. If m <= 1, an !
2364! immediate return is effected. !
2365! !
2366! N (input) INTEGER !
2367! The number of columns of the matrix A. If n <= 1, an !
2368! immediate return is effected. !
2369! !
2370! C, S (input) DOUBLE PRECISION arrays, dimension !
2371! (M-1) if SIDE = 'L' !
2372! (N-1) if SIDE = 'R' !
2373! c(k) and s(k) contain the cosine and sine that define the !
2374! matrix P(k). The two by two plane rotation part of the !
2375! matrix P(k), R(k), is assumed to be of the form !
2376! R( k ) = ( c( k ) s( k ) ). !
2377! ( -s( k ) c( k ) ) !
2378! !
2379! A (input/output) DOUBLE PRECISION array, dimension (LDA,N) !
2380! The m by n matrix A. On exit, A is overwritten by P*A if !
2381! SIDE = 'R' or by A*P' if SIDE = 'L'. !
2382! !
2383! LDA (input) INTEGER !
2384! The leading dimension of the array A. LDA >= max(1,M). !
2385! !
2386! -- LAPACK auxiliary routine (version 2.0) -- !
2387! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., !
2388! Courant Institute, Argonne National Lab, and Rice University !
2389! October 31, 1992 !
2390! !
2391!=======================================================================
2392!
2393! Imported variable declarations.
2394!
2395 integer, intent(in) :: lda, m, n
2396!
2397 real (dp), intent(in ) :: c(*), s(*)
2398 real (dp), intent(inout) :: a(lda, *)
2399!
2400 character (len=1), intent(in) :: direct, pivot, side
2401!
2402! Local variable declarations.
2403!
2404 real (r8), parameter :: zero = 0.0_dp
2405 real (r8), parameter :: one = 1.0_dp
2406!
2407 integer :: i, info, j
2408!
2409 real (dp) :: ctemp, stemp, temp
2410!
2411!-----------------------------------------------------------------------
2412! Executable Statements.
2413!-----------------------------------------------------------------------
2414!
2415! Test the input parameters.
2416!
2417 info = 0
2418 IF (.NOT. (lsame(side, 'L') .OR. lsame(side, 'R'))) THEN
2419 info = 1
2420 ELSE IF (.NOT.(lsame(pivot, 'V') .OR. lsame( pivot, 'T') .OR. &
2421 & lsame(pivot, 'B'))) THEN
2422 info = 2
2423 ELSE IF (.NOT.(lsame(direct, 'F') .OR. lsame(direct, 'B'))) THEN
2424 info = 3
2425 ELSE IF (m .LT. 0) THEN
2426 info = 4
2427 ELSE IF (n .LT. 0) THEN
2428 info = 5
2429 ELSE IF (lda .LT. max(1, m)) THEN
2430 info = 9
2431 END IF
2432 IF (info .NE. 0) THEN
2433 CALL xerbla ('DLASR ', info)
2434 RETURN
2435 END IF
2436!
2437! Quick return if possible.
2438!
2439 IF ((m.EQ.0) .OR. (n.EQ.0)) RETURN
2440!
2441 IF (lsame( side, 'L')) THEN
2442!
2443! Form P * A.
2444!
2445 IF (lsame(pivot, 'V')) THEN
2446 IF (lsame(direct, 'F')) THEN
2447 DO j = 1, m - 1
2448 ctemp = c(j)
2449 stemp = s(j)
2450 IF ((ctemp.NE.one) .OR. (stemp.NE.zero)) THEN
2451 DO i = 1, n
2452 temp = a(j+1,i)
2453 a( j+1,i) = ctemp*temp - stemp*a(j,i)
2454 a( j, i) = stemp*temp + ctemp*a(j,i)
2455 END DO
2456 END IF
2457 END DO
2458 ELSE IF (lsame(direct, 'B')) THEN
2459 DO j = m - 1, 1, -1
2460 ctemp = c(j)
2461 stemp = s(j)
2462 IF ((ctemp.NE.one) .OR. ( stemp.NE.zero ) ) THEN
2463 DO i = 1, n
2464 temp = a(j+1,i)
2465 a(j+1,i) = ctemp*temp - stemp*a(j,i)
2466 a(j, i) = stemp*temp + ctemp*a(j,i)
2467 END DO
2468 END IF
2469 END DO
2470 END IF
2471 ELSE IF (lsame(pivot, 'T')) THEN
2472 IF (lsame(direct, 'F')) THEN
2473 DO j = 2, m
2474 ctemp = c(j-1)
2475 stemp = s(j-1)
2476 IF ((ctemp.NE.one) .OR. (stemp.NE.zero)) THEN
2477 DO i = 1, n
2478 temp = a(j,i)
2479 a(j,i) = ctemp*temp - stemp*a(1,i)
2480 a(1,i) = stemp*temp + ctemp*a(1,i)
2481 END DO
2482 END IF
2483 END DO
2484 ELSE IF (lsame( direct, 'B')) THEN
2485 DO j = m, 2, -1
2486 ctemp = c(j-1)
2487 stemp = s(j-1)
2488 IF ((ctemp.NE.one) .OR. (stemp.NE.zero)) THEN
2489 DO i = 1, n
2490 temp = a(j,i)
2491 a(j,i) = ctemp*temp - stemp*a(1,i)
2492 a(1,i) = stemp*temp + ctemp*a(1,i)
2493 END DO
2494 END IF
2495 END DO
2496 END IF
2497 ELSE IF (lsame(pivot, 'B')) THEN
2498 IF (lsame(direct, 'F')) THEN
2499 DO j = 1, m - 1
2500 ctemp = c(j)
2501 stemp = s(j)
2502 IF ((ctemp.NE.one) .OR. (stemp.NE.zero)) THEN
2503 DO i = 1, n
2504 temp = a(j,i)
2505 a(j,i) = stemp*a(m,i) + ctemp*temp
2506 a(m,i) = ctemp*a(m,i) - stemp*temp
2507 END DO
2508 END IF
2509 END DO
2510 ELSE IF (lsame(direct, 'B')) THEN
2511 DO j = m - 1, 1, -1
2512 ctemp = c(j)
2513 stemp = s(j)
2514 IF ((ctemp.NE.one) .OR. (stemp.NE.zero)) THEN
2515 DO i = 1, n
2516 temp = a(j,i)
2517 a(j,i) = stemp*a(m,i) + ctemp*temp
2518 a(m,i) = ctemp*a(m,i) - stemp*temp
2519 END DO
2520 END IF
2521 END DO
2522 END IF
2523 END IF
2524 ELSE IF (lsame(side, 'R')) THEN
2525!
2526! Form A * P'.
2527!
2528 IF (lsame(pivot, 'V')) THEN
2529 IF (lsame(direct, 'F')) THEN
2530 DO j = 1, n - 1
2531 ctemp = c(j)
2532 stemp = s(j)
2533 IF ((ctemp.NE.one) .OR. (stemp.NE.zero)) THEN
2534 DO i = 1, m
2535 temp = a(i,j+1)
2536 a(i,j+1) = ctemp*temp - stemp*a(i,j)
2537 a(i,j ) = stemp*temp + ctemp*a(i,j)
2538 END DO
2539 END IF
2540 END DO
2541 ELSE IF (lsame(direct, 'B')) THEN
2542 DO j = n - 1, 1, -1
2543 ctemp = c(j)
2544 stemp = s(j)
2545 IF ((ctemp.NE.one) .OR. (stemp.NE.zero)) THEN
2546 DO i = 1, m
2547 temp = a(i,j+1)
2548 a(i,j+1) = ctemp*temp - stemp*a(i,j)
2549 a(i,j ) = stemp*temp + ctemp*a(i,j)
2550 END DO
2551 END IF
2552 END DO
2553 END IF
2554 ELSE IF (lsame(pivot, 'T')) THEN
2555 IF (lsame(direct, 'F')) THEN
2556 DO j = 2, n
2557 ctemp = c(j-1)
2558 stemp = s(j-1)
2559 IF ((ctemp.NE.one) .OR. (stemp.NE.zero)) THEN
2560 DO i = 1, m
2561 temp = a(i,j)
2562 a(i,j) = ctemp*temp - stemp*a(i,1)
2563 a(i,1) = stemp*temp + ctemp*a(i,1)
2564 END DO
2565 END IF
2566 END DO
2567 ELSE IF (lsame(direct, 'B')) THEN
2568 DO j = n, 2, -1
2569 ctemp = c(j-1)
2570 stemp = s(j-1)
2571 IF ((ctemp.NE.one) .OR. (stemp.NE.zero)) THEN
2572 DO i = 1, m
2573 temp = a(i,j)
2574 a(i,j) = ctemp*temp - stemp*a(i,1)
2575 a(i,1) = stemp*temp + ctemp*a(i,1)
2576 END DO
2577 END IF
2578 END DO
2579 END IF
2580 ELSE IF (lsame(pivot, 'B')) THEN
2581 IF (lsame(direct, 'F')) THEN
2582 DO j = 1, n - 1
2583 ctemp = c(j)
2584 stemp = s(j)
2585 IF ((ctemp.NE.one ) .OR. (stemp.NE.zero)) THEN
2586 DO i = 1, m
2587 temp = a(i,j)
2588 a(i,j) = stemp*a(i,n) + ctemp*temp
2589 a(i,n) = ctemp*a(i,n) - stemp*temp
2590 END DO
2591 END IF
2592 END DO
2593 ELSE IF (lsame(direct, 'B')) THEN
2594 DO j = n - 1, 1, -1
2595 ctemp = c(j)
2596 stemp = s(j)
2597 IF ((ctemp.NE.one) .OR. (stemp.NE.zero)) THEN
2598 DO i = 1, m
2599 temp = a(i,j)
2600 a(i,j) = stemp*a(i,n) + ctemp*temp
2601 a(i,n) = ctemp*a(i,n) - stemp*temp
2602 END DO
2603 END IF
2604 END DO
2605 END IF
2606 END IF
2607 END IF
2608!
2609 RETURN
2610 END SUBROUTINE dlasr
2611!
2612 SUBROUTINE dlasrt (ID, N, D, INFO)
2613!
2614!=======================================================================
2615! !
2616! Sort the numbers in D in increasing order (if ID = 'I') or !
2617! in decreasing order (if ID = 'D' ). !
2618! !
2619! Use Quick Sort, reverting to Insertion sort on arrays of !
2620! size <= 20. Dimension of STACK limits N to about 2**32. !
2621! !
2622! Arguments: !
2623! !
2624! ID (input) CHARACTER*1 !
2625! = 'I': sort D in increasing order; !
2626! = 'D': sort D in decreasing order. !
2627! !
2628! N (input) INTEGER !
2629! The length of the array D. !
2630! !
2631! D (input/output) DOUBLE PRECISION array, dimension (N) !
2632! On entry, the array to be sorted. !
2633! On exit, D has been sorted into increasing order !
2634! (D(1) <= ... <= D(N) ) or into decreasing order !
2635! (D(1) >= ... >= D(N) ), depending on ID. !
2636! !
2637! INFO (output) INTEGER !
2638! = 0: successful exit !
2639! < 0: if INFO = -i, the i-th argument had an illegal value !
2640! !
2641! -- LAPACK routine (version 2.0) -- !
2642! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., !
2643! Courant Institute, Argonne National Lab, and Rice University !
2644! September 30, 1994 !
2645! !
2646!=======================================================================
2647!
2648! Imported variable declarations.
2649!
2650 integer, intent(in ) :: n
2651 integer, intent(out) :: info
2652!
2653 real (dp), intent(inout) :: d(*)
2654!
2655 character (len=1), intent(in) :: id
2656!
2657! Local variable declarations.
2658!
2659 integer, parameter :: select = 20
2660!
2661 integer :: dir, endd, i, j, start, stkpnt
2662 integer :: stack( 2, 32 )
2663!
2664 real (dp) :: d1, d2, d3, dmnmx, tmp
2665!
2666!-----------------------------------------------------------------------
2667! Executable Statements.
2668!-----------------------------------------------------------------------
2669!
2670! Test the input paramters.
2671!
2672 info = 0
2673 dir = -1
2674 IF (lsame(id, 'D')) THEN
2675 dir = 0
2676 ELSE IF (lsame(id, 'I')) THEN
2677 dir = 1
2678 END IF
2679!
2680 IF (dir .EQ. -1) THEN
2681 info = -1
2682 ELSE IF (n .LT. 0) THEN
2683 info = -2
2684 END IF
2685 IF (info .NE. 0) THEN
2686 CALL xerbla ('DLASRT', -info)
2687 RETURN
2688 END IF
2689!
2690! Quick return if possible
2691!
2692 IF (n .LE. 1) RETURN
2693!
2694 stkpnt = 1
2695 stack(1,1) = 1
2696 stack(2,1) = n
2697!
2698 sort_loop : DO WHILE (stkpnt .GT. 0)
2699 start = stack(1,stkpnt)
2700 endd = stack(2,stkpnt)
2701 stkpnt = stkpnt - 1
2702 IF ((endd-start.LE.select) .AND. (endd-start.GT.0)) THEN
2703!
2704! Do Insertion sort on D(START:ENDD).
2705!
2706 IF (dir .EQ. 0) THEN
2707!
2708! Sort into decreasing order.
2709!
2710 outer1 : DO i = start + 1, endd
2711 inner1 : DO j = i, start + 1, -1
2712 IF (d(j) .GT. d(j-1)) THEN
2713 dmnmx = d(j)
2714 d(j ) = d(j-1)
2715 d(j-1) = dmnmx
2716 ELSE
2717 EXIT inner1
2718 END IF
2719 END DO inner1
2720 END DO outer1
2721!
2722 ELSE
2723!
2724! Sort into increasing order.
2725!
2726 outer2 : DO i = start + 1, endd
2727 inner2 : DO j = i, start + 1, -1
2728 IF (d(j) .LT. d(j-1)) THEN
2729 dmnmx = d(j)
2730 d(j ) = d(j-1)
2731 d(j-1) = dmnmx
2732 ELSE
2733 EXIT inner2
2734 END IF
2735 END DO inner2
2736 END DO outer2
2737!
2738 END IF
2739!
2740 ELSE IF (endd-start .GT. select) THEN
2741!
2742! Partition D( START:ENDD ) and stack parts, largest one first
2743!
2744! Choose partition entry as median of 3
2745!
2746 d1 = d(start)
2747 d2 = d(endd)
2748 i = (start+endd) / 2
2749 d3 = d( i )
2750 IF (d1 .LT. d2) THEN
2751 IF (d3 .LT. d1) THEN
2752 dmnmx = d1
2753 ELSE IF (d3 .LT. d2) THEN
2754 dmnmx = d3
2755 ELSE
2756 dmnmx = d2
2757 END IF
2758 ELSE
2759 IF (d3 .LT. d2) THEN
2760 dmnmx = d2
2761 ELSE IF (d3 .LT. d1) THEN
2762 dmnmx = d3
2763 ELSE
2764 dmnmx = d1
2765 END IF
2766 END IF
2767!
2768 IF (dir .EQ. 0) THEN
2769!
2770! Sort into decreasing order.
2771!
2772 i = start - 1
2773 j = endd + 1
2774#ifdef REMOVE_LAPACK_GOTOS
2775 decloop_3 : DO ! forever LOOP 3
2776 decloop_1 : DO ! forever LOOP 1
2777 j = j - 1
2778 IF (d(j) .GE. dmnmx) EXIT decloop_1 ! terminate LOOP 1
2779 END DO decloop_1
2780!
2781 decloop_2 : DO ! forever LOOP 2
2782 i = i + 1
2783 IF (d(i) .LE. dmnmx) EXIT decloop_2 ! terminate LOOP 2
2784 END DO decloop_2
2785!
2786 IF (i .LT. j) THEN
2787 tmp = d(i)
2788 d(i) = d(j)
2789 d(j) = tmp
2790 cycle decloop_3 ! iterate LOOP 3
2791 ELSE
2792 EXIT decloop_3 ! terminate LOOP 3
2793 END IF
2794 END DO decloop_3
2795#else
2796 60 CONTINUE
2797 70 CONTINUE
2798 j = j - 1
2799 IF (d(j) .LT. dmnmx) GO TO 70
2800 80 CONTINUE
2801 i = i + 1
2802 IF (d(i) .GT. dmnmx) GO TO 80
2803 IF (i .LT. j) THEN
2804 tmp = d(i)
2805 d(i) = d(j)
2806 d(j) = tmp
2807 GO TO 60
2808 END IF
2809#endif
2810 IF (j-start .GT. endd-j-1) THEN
2811 stkpnt = stkpnt + 1
2812 stack(1, stkpnt) = start
2813 stack(2, stkpnt) = j
2814 stkpnt = stkpnt + 1
2815 stack(1, stkpnt) = j + 1
2816 stack(2, stkpnt) = endd
2817 ELSE
2818 stkpnt = stkpnt + 1
2819 stack(1, stkpnt) = j + 1
2820 stack(2, stkpnt) = endd
2821 stkpnt = stkpnt + 1
2822 stack(1, stkpnt) = start
2823 stack(2, stkpnt) = j
2824 END IF
2825 ELSE
2826!
2827! Sort into increasing order.
2828!
2829 i = start - 1
2830 j = endd + 1
2831#ifdef REMOVE_LAPACK_GOTOS
2832 incloop_3 : DO ! forever LOOP 3
2833 incloop_1 : DO ! forever LOOP 1
2834 j = j - 1
2835 IF (d(j) .LE. dmnmx) EXIT incloop_1 ! terminate LOOP 1
2836 END DO incloop_1
2837!
2838 incloop_2 : DO ! forever LOOP 2
2839 i = i + 1
2840 IF (d(i) .GE. dmnmx) EXIT incloop_2 ! terminate LOOP 2
2841 END DO incloop_2
2842!
2843 IF (i .LT. j) THEN
2844 tmp = d(i)
2845 d(i) = d(j)
2846 d(j) = tmp
2847 cycle incloop_3 ! iterate LOOP 3
2848 ELSE
2849 EXIT incloop_3 ! terminate LOOP 3
2850 END IF
2851 END DO incloop_3
2852#else
2853 90 CONTINUE
2854 100 CONTINUE
2855 j = j - 1
2856 IF (d(j) .GT. dmnmx) GO TO 100
2857 110 CONTINUE
2858 i = i + 1
2859 IF (d(i) .LT. dmnmx) GO TO 110
2860 IF (i .LT. j) THEN
2861 tmp = d(i)
2862 d(i) = d(j)
2863 d(j) = tmp
2864 GO TO 90
2865 END IF
2866#endif
2867 IF (j-start .GT. endd-j-1) THEN
2868 stkpnt = stkpnt + 1
2869 stack(1, stkpnt) = start
2870 stack(2, stkpnt) = j
2871 stkpnt = stkpnt + 1
2872 stack(1, stkpnt) = j + 1
2873 stack(2, stkpnt) = endd
2874 ELSE
2875 stkpnt = stkpnt + 1
2876 stack(1, stkpnt) = j + 1
2877 stack(2, stkpnt) = endd
2878 stkpnt = stkpnt + 1
2879 stack(1, stkpnt) = start
2880 stack(2, stkpnt) = j
2881 END IF
2882 END IF
2883 END IF
2884 END DO sort_loop
2885!
2886 RETURN
2887 END SUBROUTINE dlasrt
2888!
2889 SUBROUTINE dlassq (N, X, INCX, SCALE, SUMSQ)
2890!
2891!=======================================================================
2892! !
2893! DLASSQ returns the values scl and smsq such that !
2894! !
2895! ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq !
2896! !
2897! where x( i ) = X( 1 + ( i - 1 )*INCX ). The value of sumsq is !
2898! assumed to be non-negative and scl returns the value !
2899! !
2900! scl = max( scale, abs( x( i ) ) ). !
2901! !
2902! scale and sumsq must be supplied in SCALE and SUMSQ and !
2903! scl and smsq are overwritten on SCALE and SUMSQ respectively. !
2904! !
2905! The routine makes only one pass through the vector x. !
2906! !
2907! Arguments: !
2908! !
2909! N (input) INTEGER !
2910! The number of elements to be used from the vector X. !
2911! !
2912! X (input) DOUBLE PRECISION !
2913! The vector for which a scaled sum of squares is computed. !
2914! x( i ) = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n. !
2915! !
2916! INCX (input) INTEGER !
2917! The increment between successive values of the vector X. !
2918! INCX > 0. !
2919! !
2920! SCALE (input/output) DOUBLE PRECISION !
2921! On entry, the value scale in the equation above. !
2922! On exit, SCALE is overwritten with scl, the scaling factor !
2923! for the sum of squares. !
2924! !
2925! SUMSQ (input/output) DOUBLE PRECISION !
2926! On entry, the value sumsq in the equation above. !
2927! On exit, SUMSQ is overwritten with smsq , the basic sum of !
2928! squares from which scl has been factored out. !
2929! !
2930! -- LAPACK auxiliary routine (version 2.0) -- !
2931! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., !
2932! Courant Institute, Argonne National Lab, and Rice University !
2933! October 31, 1992 !
2934! !
2935!=======================================================================
2936!
2937! Imported variable declarations.
2938!
2939 integer, intent(in) :: incx, n
2940!
2941 real (dp), intent(in ) :: x(*)
2942 real (dp), intent(inout) :: scale, sumsq
2943!
2944! Local variable declarations.
2945!
2946 real (dp), parameter :: zero = 0.0_dp
2947!
2948 integer :: ix
2949!
2950 real (dp) :: absxi
2951!
2952!-----------------------------------------------------------------------
2953! Executable Statements.
2954!-----------------------------------------------------------------------
2955!
2956 IF (n .GT. 0) THEN
2957 DO ix = 1, 1 + ( n-1 )*incx, incx
2958 IF (x(ix) .NE. zero) THEN
2959 absxi = abs(x(ix))
2960 IF (scale .LT. absxi) THEN
2961 sumsq = 1.0_dp + sumsq*(scale / absxi)**2
2962 scale = absxi
2963 ELSE
2964 sumsq = sumsq + (absxi / scale)**2
2965 END IF
2966 END IF
2967 END DO
2968 END IF
2969!
2970 RETURN
2971 END SUBROUTINE dlassq
2972!
2973 SUBROUTINE dswap (N, DX, INCX, DY, INCY)
2974!
2975!=======================================================================
2976! !
2977! DSWAP interchanges two vectors. It uses unrolled loops for !
2978! increments equal one. !
2979! !
2980! Arguments: !
2981! !
2982! N (input) INTEGER !
2983! number of elements in input vector(s) !
2984! !
2985! DX (input/output) DOUBLE PRECISION !
2986! array, dimension ( 1 + ( N - 1 )*abs( INCX ) ) !
2987! !
2988! INCX (input) INTEGER !
2989! storage spacing between elements of DX !
2990! !
2991! DY (input/output) DOUBLE PRECISION !
2992! array, dimension ( 1 + ( N - 1 )*abs( INCY ) ) !
2993! !
2994! INCY (input) INTEGER !
2995! storage spacing between elements of DY !
2996! !
2997! -- LAPACK auxiliary routine (version 2.0) -- !
2998! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., !
2999! Courant Institute, Argonne National Lab, and Rice University !
3000! October 31, 1992 !
3001! !
3002! Based on Jack Dongarra, LINPACK, 3/11/78. !
3003! !
3004!=======================================================================
3005!
3006! Imported variable declarations.
3007!
3008 integer, intent(in) :: INCX, INCY, N
3009!
3010 real (dp), intent(inout) :: DX(*), DY(*)
3011!
3012! Local variable declarations.
3013!
3014 integer I, IX, IY, M, MP1
3015!
3016 real (dp) :: DTEMP
3017!
3018!-----------------------------------------------------------------------
3019! Executable Statements.
3020!-----------------------------------------------------------------------
3021!
3022 IF (n .le. 0) RETURN
3023!
3024 IF (incx.EQ.1 .AND. incy.EQ.1) THEN
3025!
3026! Code for both increments equal to 1.
3027!
3028 m = mod(n,3)
3029 IF (m .NE. 0 ) THEN
3030 DO i = 1, m
3031 dtemp = dx(i)
3032 dx(i) = dy(i)
3033 dy(i) = dtemp
3034 END DO
3035 IF (n .LT. 3) RETURN
3036 END IF
3037 mp1 = m + 1
3038 DO i = mp1, n, 3
3039 dtemp = dx(i)
3040 dx(i) = dy(i)
3041 dy(i) = dtemp
3042 dtemp = dx(i+1)
3043 dx(i+1) = dy(i+1)
3044 dy(i+1) = dtemp
3045 dtemp = dx(i+2)
3046 dx(i+2) = dy(i+2)
3047 dy(i+2) = dtemp
3048 END DO
3049 ELSE
3050!
3051! Code for unequal increments or equal increments not equal to 1.
3052!
3053 ix = 1
3054 iy = 1
3055 IF (incx .LT. 0) ix = (-n+1)*incx + 1
3056 IF (incy .LT. 0) iy = (-n+1)*incy + 1
3057 DO i = 1, n
3058 dtemp = dx(ix)
3059 dx(ix) = dy(iy)
3060 dy(iy) = dtemp
3061 ix = ix + incx
3062 iy = iy + incy
3063 END DO
3064 END IF
3065!
3066 RETURN
3067 END SUBROUTINE dswap
3068!
3069 FUNCTION lsame (CA, CB) RESULT (IsSAME)
3070!
3071!=======================================================================
3072! !
3073! LSAME returns .TRUE. if CA is the same letter as CB regardless of !
3074! case. !
3075! !
3076! Arguments: !
3077! !
3078! CA (input) CHARACTER*1 !
3079! CB (input) CHARACTER*1 !
3080! CA and CB specify the single characters to be compared. !
3081! !
3082! -- LAPACK auxiliary routine (version 2.0) -- !
3083! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., !
3084! Courant Institute, Argonne National Lab, and Rice University !
3085! September 30, 1994 !
3086! !
3087!=======================================================================
3088!
3089! Imported variable declarations.
3090!
3091 character (len=1), intent(in) :: ca, cb
3092!
3093! Local variable declarations.
3094!
3095 logical :: issame
3096!
3097 integer :: inta, intb, zcode
3098!
3099!-----------------------------------------------------------------------
3100! Executable Statements.
3101!-----------------------------------------------------------------------
3102!
3103! Test if the characters are equal.
3104!
3105 issame = ca .EQ. cb
3106 IF( issame ) RETURN
3107!
3108! Now test for equivalence if both characters are alphabetic.
3109!
3110 zcode = ichar('Z')
3111!
3112! Use 'Z' rather than 'A' so that ASCII can be detected on Prime
3113! machines, on which ICHAR returns a value with bit 8 set.
3114! ICHAR('A') on Prime machines returns 193 which is the same as
3115! ICHAR('A') on an EBCDIC machine.
3116!
3117 inta = ichar(ca)
3118 intb = ichar(cb)
3119!
3120 IF (zcode.EQ.90 .OR. zcode.EQ.122) THEN
3121!
3122! ASCII is assumed - ZCODE is the ASCII code of either lower or
3123! upper case 'Z'.
3124!
3125 IF (inta.GE.97 .AND. inta.LE.122) inta = inta - 32
3126 IF (intb.GE.97 .AND. intb.LE.122) intb = intb - 32
3127!
3128 ELSE IF (zcode.EQ.233 .OR. zcode.EQ.169) THEN
3129!
3130! EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or
3131! upper case 'Z'.
3132!
3133 IF (inta.GE.129 .AND. inta.LE.137 .OR. &
3134 & inta.GE.145 .AND. inta.LE.153 .OR. &
3135 & inta.GE.162 .AND. inta.LE.169) inta = inta + 64
3136 IF (intb.GE.129 .AND. intb.LE.137 .OR. &
3137 & intb.GE.145 .AND. intb.LE.153 .OR. &
3138 & intb.GE.162 .AND. intb.LE.169) intb = intb + 64
3139!
3140 ELSE IF (zcode.EQ.218 .OR. zcode.EQ.250) THEN
3141!
3142! ASCII is assumed, on Prime machines - ZCODE is the ASCII code
3143! plus 128 of either lower or upper case 'Z'.
3144!
3145 IF (inta.GE.225 .AND. inta.LE.250) inta = inta - 32
3146 IF (intb.GE.225 .AND. intb.LE.250) intb = intb - 32
3147 END IF
3148!
3149 issame = inta .EQ. intb
3150
3151 RETURN
3152 END FUNCTION lsame
3153!
3154 SUBROUTINE xerbla (SRNAME, INFO)
3155!
3156!=======================================================================
3157! !
3158! XERBLA is an error handler for the LAPACK routines. !
3159! It is called by an LAPACK routine if an input parameter has an !
3160! invalid value. A message is printed and execution stops. !
3161! !
3162! Installers may consider modifying the STOP statement in order to !
3163! call system-specific exception-handling facilities. !
3164! !
3165! Arguments: !
3166! !
3167! SRNAME (input) CHARACTER*6 !
3168! The name of the routine which called XERBLA. !
3169! !
3170! INFO (input) INTEGER !
3171! The position of the invalid parameter in the parameter list !
3172! of the calling routine. !
3173! !
3174! -- LAPACK auxiliary routine (version 2.0) -- !
3175! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., !
3176! Courant Institute, Argonne National Lab, and Rice University !
3177! September 30, 1994 !
3178! !
3179!=======================================================================
3180!
3181! Imported variable declarations.
3182!
3183 integer, intent(in) :: info
3184!
3185 character (len=*), intent(in) :: srname
3186!
3187!-----------------------------------------------------------------------
3188! Executable Statements.
3189!-----------------------------------------------------------------------
3190!
3191 print 10, srname, info
3192 10 FORMAT (' LAPACK_MOD: ** On entry to ', a, ' parameter number ', &
3193 & i0,' had an illegal value.')
3194 stop
3195!
3196 END SUBROUTINE xerbla
3197!
3198 END MODULE lapack_mod
lapack_mod
Definition lapack_mod.F:2
lapack_mod::dlanst
function, private dlanst(norm, n, d, e)
Definition lapack_mod.F:1605
lapack_mod::dlaev2
subroutine, private dlaev2(a, b, c, rt1, rt2, cs1, sn1)
Definition lapack_mod.F:651
lapack_mod::dlasr
subroutine, private dlasr(side, pivot, direct, m, n, c, s, a, lda)
Definition lapack_mod.F:2298
lapack_mod::dlamc5
subroutine, private dlamc5(beta, p, emin, ieee, emax, rmax)
Definition lapack_mod.F:1325
lapack_mod::dlamc3
real(dp) function, private dlamc3(a, b)
Definition lapack_mod.F:1213
lapack_mod::dswap
subroutine dswap(n, dx, incx, dy, incy)
Definition lapack_mod.F:2974
lapack_mod::dlapy2
real(dp) function, private dlapy2(x, y)
Definition lapack_mod.F:1719
lapack_mod::dlamc1
subroutine, private dlamc1(beta, t, rnd, ieee1)
Definition lapack_mod.F:813
lapack_mod::dlartg
subroutine, private dlartg(f, g, cs, sn, r)
Definition lapack_mod.F:1768
lapack_mod::dlasrt
subroutine, private dlasrt(id, n, d, info)
Definition lapack_mod.F:2613
lapack_mod::xerbla
subroutine, private xerbla(srname, info)
Definition lapack_mod.F:3155
lapack_mod::dlamch
real(dp) function, private dlamch(cmach)
Definition lapack_mod.F:1491
lapack_mod::dlascl
subroutine, private dlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
Definition lapack_mod.F:1931
lapack_mod::lsame
logical function, private lsame(ca, cb)
Definition lapack_mod.F:3070
lapack_mod::dlamc4
subroutine, private dlamc4(emin, start, base)
Definition lapack_mod.F:1250
lapack_mod::dlaset
subroutine, private dlaset(uplo, m, n, alpha, beta, a, lda)
Definition lapack_mod.F:2188
lapack_mod::dlae2
subroutine, private dlae2(a, b, c, rt1, rt2)
Definition lapack_mod.F:534
lapack_mod::dlassq
subroutine, private dlassq(n, x, incx, scale, sumsq)
Definition lapack_mod.F:2890
lapack_mod::dsteqr
subroutine, public dsteqr(compz, n, d, e, z, ldz, work, info)
Definition lapack_mod.F:66
lapack_mod::dlamc2
subroutine, private dlamc2(beta, t, rnd, eps, emin, rmin, emax, rmax)
Definition lapack_mod.F:976
mod_kinds
Definition mod_kinds.F:2
mod_kinds::r8
integer, parameter r8
Definition mod_kinds.F:28
mod_kinds::dp
integer, parameter dp
Definition mod_kinds.F:25