source: ThirdParty/mpqc_open/src/lib/math/scmat/pdsteqr.f@ e9f307

Candidate_v1.7.0 stable
Last change on this file since e9f307 was 860145, checked in by Frederik Heber <heber@…>, 9 years ago

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