1: SUBROUTINE DLALN2( LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B,
2: $ LDB, WR, WI, X, LDX, SCALE, XNORM, INFO )
3: *
4: * -- LAPACK auxiliary routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
8: *
9: * .. Scalar Arguments ..
10: LOGICAL LTRANS
11: INTEGER INFO, LDA, LDB, LDX, NA, NW
12: DOUBLE PRECISION CA, D1, D2, SCALE, SMIN, WI, WR, XNORM
13: * ..
14: * .. Array Arguments ..
15: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), X( LDX, * )
16: * ..
17: *
18: * Purpose
19: * =======
20: *
21: * DLALN2 solves a system of the form (ca A - w D ) X = s B
22: * or (ca A' - w D) X = s B with possible scaling ("s") and
23: * perturbation of A. (A' means A-transpose.)
24: *
25: * A is an NA x NA real matrix, ca is a real scalar, D is an NA x NA
26: * real diagonal matrix, w is a real or complex value, and X and B are
27: * NA x 1 matrices -- real if w is real, complex if w is complex. NA
28: * may be 1 or 2.
29: *
30: * If w is complex, X and B are represented as NA x 2 matrices,
31: * the first column of each being the real part and the second
32: * being the imaginary part.
33: *
34: * "s" is a scaling factor (.LE. 1), computed by DLALN2, which is
35: * so chosen that X can be computed without overflow. X is further
36: * scaled if necessary to assure that norm(ca A - w D)*norm(X) is less
37: * than overflow.
38: *
39: * If both singular values of (ca A - w D) are less than SMIN,
40: * SMIN*identity will be used instead of (ca A - w D). If only one
41: * singular value is less than SMIN, one element of (ca A - w D) will be
42: * perturbed enough to make the smallest singular value roughly SMIN.
43: * If both singular values are at least SMIN, (ca A - w D) will not be
44: * perturbed. In any case, the perturbation will be at most some small
45: * multiple of max( SMIN, ulp*norm(ca A - w D) ). The singular values
46: * are computed by infinity-norm approximations, and thus will only be
47: * correct to a factor of 2 or so.
48: *
49: * Note: all input quantities are assumed to be smaller than overflow
50: * by a reasonable factor. (See BIGNUM.)
51: *
52: * Arguments
53: * ==========
54: *
55: * LTRANS (input) LOGICAL
56: * =.TRUE.: A-transpose will be used.
57: * =.FALSE.: A will be used (not transposed.)
58: *
59: * NA (input) INTEGER
60: * The size of the matrix A. It may (only) be 1 or 2.
61: *
62: * NW (input) INTEGER
63: * 1 if "w" is real, 2 if "w" is complex. It may only be 1
64: * or 2.
65: *
66: * SMIN (input) DOUBLE PRECISION
67: * The desired lower bound on the singular values of A. This
68: * should be a safe distance away from underflow or overflow,
69: * say, between (underflow/machine precision) and (machine
70: * precision * overflow ). (See BIGNUM and ULP.)
71: *
72: * CA (input) DOUBLE PRECISION
73: * The coefficient c, which A is multiplied by.
74: *
75: * A (input) DOUBLE PRECISION array, dimension (LDA,NA)
76: * The NA x NA matrix A.
77: *
78: * LDA (input) INTEGER
79: * The leading dimension of A. It must be at least NA.
80: *
81: * D1 (input) DOUBLE PRECISION
82: * The 1,1 element in the diagonal matrix D.
83: *
84: * D2 (input) DOUBLE PRECISION
85: * The 2,2 element in the diagonal matrix D. Not used if NW=1.
86: *
87: * B (input) DOUBLE PRECISION array, dimension (LDB,NW)
88: * The NA x NW matrix B (right-hand side). If NW=2 ("w" is
89: * complex), column 1 contains the real part of B and column 2
90: * contains the imaginary part.
91: *
92: * LDB (input) INTEGER
93: * The leading dimension of B. It must be at least NA.
94: *
95: * WR (input) DOUBLE PRECISION
96: * The real part of the scalar "w".
97: *
98: * WI (input) DOUBLE PRECISION
99: * The imaginary part of the scalar "w". Not used if NW=1.
100: *
101: * X (output) DOUBLE PRECISION array, dimension (LDX,NW)
102: * The NA x NW matrix X (unknowns), as computed by DLALN2.
103: * If NW=2 ("w" is complex), on exit, column 1 will contain
104: * the real part of X and column 2 will contain the imaginary
105: * part.
106: *
107: * LDX (input) INTEGER
108: * The leading dimension of X. It must be at least NA.
109: *
110: * SCALE (output) DOUBLE PRECISION
111: * The scale factor that B must be multiplied by to insure
112: * that overflow does not occur when computing X. Thus,
113: * (ca A - w D) X will be SCALE*B, not B (ignoring
114: * perturbations of A.) It will be at most 1.
115: *
116: * XNORM (output) DOUBLE PRECISION
117: * The infinity-norm of X, when X is regarded as an NA x NW
118: * real matrix.
119: *
120: * INFO (output) INTEGER
121: * An error flag. It will be set to zero if no error occurs,
122: * a negative number if an argument is in error, or a positive
123: * number if ca A - w D had to be perturbed.
124: * The possible values are:
125: * = 0: No error occurred, and (ca A - w D) did not have to be
126: * perturbed.
127: * = 1: (ca A - w D) had to be perturbed to make its smallest
128: * (or only) singular value greater than SMIN.
129: * NOTE: In the interests of speed, this routine does not
130: * check the inputs for errors.
131: *
132: * =====================================================================
133: *
134: * .. Parameters ..
135: DOUBLE PRECISION ZERO, ONE
136: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
137: DOUBLE PRECISION TWO
138: PARAMETER ( TWO = 2.0D0 )
139: * ..
140: * .. Local Scalars ..
141: INTEGER ICMAX, J
142: DOUBLE PRECISION BBND, BI1, BI2, BIGNUM, BNORM, BR1, BR2, CI21,
143: $ CI22, CMAX, CNORM, CR21, CR22, CSI, CSR, LI21,
144: $ LR21, SMINI, SMLNUM, TEMP, U22ABS, UI11, UI11R,
145: $ UI12, UI12S, UI22, UR11, UR11R, UR12, UR12S,
146: $ UR22, XI1, XI2, XR1, XR2
147: * ..
148: * .. Local Arrays ..
149: LOGICAL RSWAP( 4 ), ZSWAP( 4 )
150: INTEGER IPIVOT( 4, 4 )
151: DOUBLE PRECISION CI( 2, 2 ), CIV( 4 ), CR( 2, 2 ), CRV( 4 )
152: * ..
153: * .. External Functions ..
154: DOUBLE PRECISION DLAMCH
155: EXTERNAL DLAMCH
156: * ..
157: * .. External Subroutines ..
158: EXTERNAL DLADIV
159: * ..
160: * .. Intrinsic Functions ..
161: INTRINSIC ABS, MAX
162: * ..
163: * .. Equivalences ..
164: EQUIVALENCE ( CI( 1, 1 ), CIV( 1 ) ),
165: $ ( CR( 1, 1 ), CRV( 1 ) )
166: * ..
167: * .. Data statements ..
168: DATA ZSWAP / .FALSE., .FALSE., .TRUE., .TRUE. /
169: DATA RSWAP / .FALSE., .TRUE., .FALSE., .TRUE. /
170: DATA IPIVOT / 1, 2, 3, 4, 2, 1, 4, 3, 3, 4, 1, 2, 4,
171: $ 3, 2, 1 /
172: * ..
173: * .. Executable Statements ..
174: *
175: * Compute BIGNUM
176: *
177: SMLNUM = TWO*DLAMCH( 'Safe minimum' )
178: BIGNUM = ONE / SMLNUM
179: SMINI = MAX( SMIN, SMLNUM )
180: *
181: * Don't check for input errors
182: *
183: INFO = 0
184: *
185: * Standard Initializations
186: *
187: SCALE = ONE
188: *
189: IF( NA.EQ.1 ) THEN
190: *
191: * 1 x 1 (i.e., scalar) system C X = B
192: *
193: IF( NW.EQ.1 ) THEN
194: *
195: * Real 1x1 system.
196: *
197: * C = ca A - w D
198: *
199: CSR = CA*A( 1, 1 ) - WR*D1
200: CNORM = ABS( CSR )
201: *
202: * If | C | < SMINI, use C = SMINI
203: *
204: IF( CNORM.LT.SMINI ) THEN
205: CSR = SMINI
206: CNORM = SMINI
207: INFO = 1
208: END IF
209: *
210: * Check scaling for X = B / C
211: *
212: BNORM = ABS( B( 1, 1 ) )
213: IF( CNORM.LT.ONE .AND. BNORM.GT.ONE ) THEN
214: IF( BNORM.GT.BIGNUM*CNORM )
215: $ SCALE = ONE / BNORM
216: END IF
217: *
218: * Compute X
219: *
220: X( 1, 1 ) = ( B( 1, 1 )*SCALE ) / CSR
221: XNORM = ABS( X( 1, 1 ) )
222: ELSE
223: *
224: * Complex 1x1 system (w is complex)
225: *
226: * C = ca A - w D
227: *
228: CSR = CA*A( 1, 1 ) - WR*D1
229: CSI = -WI*D1
230: CNORM = ABS( CSR ) + ABS( CSI )
231: *
232: * If | C | < SMINI, use C = SMINI
233: *
234: IF( CNORM.LT.SMINI ) THEN
235: CSR = SMINI
236: CSI = ZERO
237: CNORM = SMINI
238: INFO = 1
239: END IF
240: *
241: * Check scaling for X = B / C
242: *
243: BNORM = ABS( B( 1, 1 ) ) + ABS( B( 1, 2 ) )
244: IF( CNORM.LT.ONE .AND. BNORM.GT.ONE ) THEN
245: IF( BNORM.GT.BIGNUM*CNORM )
246: $ SCALE = ONE / BNORM
247: END IF
248: *
249: * Compute X
250: *
251: CALL DLADIV( SCALE*B( 1, 1 ), SCALE*B( 1, 2 ), CSR, CSI,
252: $ X( 1, 1 ), X( 1, 2 ) )
253: XNORM = ABS( X( 1, 1 ) ) + ABS( X( 1, 2 ) )
254: END IF
255: *
256: ELSE
257: *
258: * 2x2 System
259: *
260: * Compute the real part of C = ca A - w D (or ca A' - w D )
261: *
262: CR( 1, 1 ) = CA*A( 1, 1 ) - WR*D1
263: CR( 2, 2 ) = CA*A( 2, 2 ) - WR*D2
264: IF( LTRANS ) THEN
265: CR( 1, 2 ) = CA*A( 2, 1 )
266: CR( 2, 1 ) = CA*A( 1, 2 )
267: ELSE
268: CR( 2, 1 ) = CA*A( 2, 1 )
269: CR( 1, 2 ) = CA*A( 1, 2 )
270: END IF
271: *
272: IF( NW.EQ.1 ) THEN
273: *
274: * Real 2x2 system (w is real)
275: *
276: * Find the largest element in C
277: *
278: CMAX = ZERO
279: ICMAX = 0
280: *
281: DO 10 J = 1, 4
282: IF( ABS( CRV( J ) ).GT.CMAX ) THEN
283: CMAX = ABS( CRV( J ) )
284: ICMAX = J
285: END IF
286: 10 CONTINUE
287: *
288: * If norm(C) < SMINI, use SMINI*identity.
289: *
290: IF( CMAX.LT.SMINI ) THEN
291: BNORM = MAX( ABS( B( 1, 1 ) ), ABS( B( 2, 1 ) ) )
292: IF( SMINI.LT.ONE .AND. BNORM.GT.ONE ) THEN
293: IF( BNORM.GT.BIGNUM*SMINI )
294: $ SCALE = ONE / BNORM
295: END IF
296: TEMP = SCALE / SMINI
297: X( 1, 1 ) = TEMP*B( 1, 1 )
298: X( 2, 1 ) = TEMP*B( 2, 1 )
299: XNORM = TEMP*BNORM
300: INFO = 1
301: RETURN
302: END IF
303: *
304: * Gaussian elimination with complete pivoting.
305: *
306: UR11 = CRV( ICMAX )
307: CR21 = CRV( IPIVOT( 2, ICMAX ) )
308: UR12 = CRV( IPIVOT( 3, ICMAX ) )
309: CR22 = CRV( IPIVOT( 4, ICMAX ) )
310: UR11R = ONE / UR11
311: LR21 = UR11R*CR21
312: UR22 = CR22 - UR12*LR21
313: *
314: * If smaller pivot < SMINI, use SMINI
315: *
316: IF( ABS( UR22 ).LT.SMINI ) THEN
317: UR22 = SMINI
318: INFO = 1
319: END IF
320: IF( RSWAP( ICMAX ) ) THEN
321: BR1 = B( 2, 1 )
322: BR2 = B( 1, 1 )
323: ELSE
324: BR1 = B( 1, 1 )
325: BR2 = B( 2, 1 )
326: END IF
327: BR2 = BR2 - LR21*BR1
328: BBND = MAX( ABS( BR1*( UR22*UR11R ) ), ABS( BR2 ) )
329: IF( BBND.GT.ONE .AND. ABS( UR22 ).LT.ONE ) THEN
330: IF( BBND.GE.BIGNUM*ABS( UR22 ) )
331: $ SCALE = ONE / BBND
332: END IF
333: *
334: XR2 = ( BR2*SCALE ) / UR22
335: XR1 = ( SCALE*BR1 )*UR11R - XR2*( UR11R*UR12 )
336: IF( ZSWAP( ICMAX ) ) THEN
337: X( 1, 1 ) = XR2
338: X( 2, 1 ) = XR1
339: ELSE
340: X( 1, 1 ) = XR1
341: X( 2, 1 ) = XR2
342: END IF
343: XNORM = MAX( ABS( XR1 ), ABS( XR2 ) )
344: *
345: * Further scaling if norm(A) norm(X) > overflow
346: *
347: IF( XNORM.GT.ONE .AND. CMAX.GT.ONE ) THEN
348: IF( XNORM.GT.BIGNUM / CMAX ) THEN
349: TEMP = CMAX / BIGNUM
350: X( 1, 1 ) = TEMP*X( 1, 1 )
351: X( 2, 1 ) = TEMP*X( 2, 1 )
352: XNORM = TEMP*XNORM
353: SCALE = TEMP*SCALE
354: END IF
355: END IF
356: ELSE
357: *
358: * Complex 2x2 system (w is complex)
359: *
360: * Find the largest element in C
361: *
362: CI( 1, 1 ) = -WI*D1
363: CI( 2, 1 ) = ZERO
364: CI( 1, 2 ) = ZERO
365: CI( 2, 2 ) = -WI*D2
366: CMAX = ZERO
367: ICMAX = 0
368: *
369: DO 20 J = 1, 4
370: IF( ABS( CRV( J ) )+ABS( CIV( J ) ).GT.CMAX ) THEN
371: CMAX = ABS( CRV( J ) ) + ABS( CIV( J ) )
372: ICMAX = J
373: END IF
374: 20 CONTINUE
375: *
376: * If norm(C) < SMINI, use SMINI*identity.
377: *
378: IF( CMAX.LT.SMINI ) THEN
379: BNORM = MAX( ABS( B( 1, 1 ) )+ABS( B( 1, 2 ) ),
380: $ ABS( B( 2, 1 ) )+ABS( B( 2, 2 ) ) )
381: IF( SMINI.LT.ONE .AND. BNORM.GT.ONE ) THEN
382: IF( BNORM.GT.BIGNUM*SMINI )
383: $ SCALE = ONE / BNORM
384: END IF
385: TEMP = SCALE / SMINI
386: X( 1, 1 ) = TEMP*B( 1, 1 )
387: X( 2, 1 ) = TEMP*B( 2, 1 )
388: X( 1, 2 ) = TEMP*B( 1, 2 )
389: X( 2, 2 ) = TEMP*B( 2, 2 )
390: XNORM = TEMP*BNORM
391: INFO = 1
392: RETURN
393: END IF
394: *
395: * Gaussian elimination with complete pivoting.
396: *
397: UR11 = CRV( ICMAX )
398: UI11 = CIV( ICMAX )
399: CR21 = CRV( IPIVOT( 2, ICMAX ) )
400: CI21 = CIV( IPIVOT( 2, ICMAX ) )
401: UR12 = CRV( IPIVOT( 3, ICMAX ) )
402: UI12 = CIV( IPIVOT( 3, ICMAX ) )
403: CR22 = CRV( IPIVOT( 4, ICMAX ) )
404: CI22 = CIV( IPIVOT( 4, ICMAX ) )
405: IF( ICMAX.EQ.1 .OR. ICMAX.EQ.4 ) THEN
406: *
407: * Code when off-diagonals of pivoted C are real
408: *
409: IF( ABS( UR11 ).GT.ABS( UI11 ) ) THEN
410: TEMP = UI11 / UR11
411: UR11R = ONE / ( UR11*( ONE+TEMP**2 ) )
412: UI11R = -TEMP*UR11R
413: ELSE
414: TEMP = UR11 / UI11
415: UI11R = -ONE / ( UI11*( ONE+TEMP**2 ) )
416: UR11R = -TEMP*UI11R
417: END IF
418: LR21 = CR21*UR11R
419: LI21 = CR21*UI11R
420: UR12S = UR12*UR11R
421: UI12S = UR12*UI11R
422: UR22 = CR22 - UR12*LR21
423: UI22 = CI22 - UR12*LI21
424: ELSE
425: *
426: * Code when diagonals of pivoted C are real
427: *
428: UR11R = ONE / UR11
429: UI11R = ZERO
430: LR21 = CR21*UR11R
431: LI21 = CI21*UR11R
432: UR12S = UR12*UR11R
433: UI12S = UI12*UR11R
434: UR22 = CR22 - UR12*LR21 + UI12*LI21
435: UI22 = -UR12*LI21 - UI12*LR21
436: END IF
437: U22ABS = ABS( UR22 ) + ABS( UI22 )
438: *
439: * If smaller pivot < SMINI, use SMINI
440: *
441: IF( U22ABS.LT.SMINI ) THEN
442: UR22 = SMINI
443: UI22 = ZERO
444: INFO = 1
445: END IF
446: IF( RSWAP( ICMAX ) ) THEN
447: BR2 = B( 1, 1 )
448: BR1 = B( 2, 1 )
449: BI2 = B( 1, 2 )
450: BI1 = B( 2, 2 )
451: ELSE
452: BR1 = B( 1, 1 )
453: BR2 = B( 2, 1 )
454: BI1 = B( 1, 2 )
455: BI2 = B( 2, 2 )
456: END IF
457: BR2 = BR2 - LR21*BR1 + LI21*BI1
458: BI2 = BI2 - LI21*BR1 - LR21*BI1
459: BBND = MAX( ( ABS( BR1 )+ABS( BI1 ) )*
460: $ ( U22ABS*( ABS( UR11R )+ABS( UI11R ) ) ),
461: $ ABS( BR2 )+ABS( BI2 ) )
462: IF( BBND.GT.ONE .AND. U22ABS.LT.ONE ) THEN
463: IF( BBND.GE.BIGNUM*U22ABS ) THEN
464: SCALE = ONE / BBND
465: BR1 = SCALE*BR1
466: BI1 = SCALE*BI1
467: BR2 = SCALE*BR2
468: BI2 = SCALE*BI2
469: END IF
470: END IF
471: *
472: CALL DLADIV( BR2, BI2, UR22, UI22, XR2, XI2 )
473: XR1 = UR11R*BR1 - UI11R*BI1 - UR12S*XR2 + UI12S*XI2
474: XI1 = UI11R*BR1 + UR11R*BI1 - UI12S*XR2 - UR12S*XI2
475: IF( ZSWAP( ICMAX ) ) THEN
476: X( 1, 1 ) = XR2
477: X( 2, 1 ) = XR1
478: X( 1, 2 ) = XI2
479: X( 2, 2 ) = XI1
480: ELSE
481: X( 1, 1 ) = XR1
482: X( 2, 1 ) = XR2
483: X( 1, 2 ) = XI1
484: X( 2, 2 ) = XI2
485: END IF
486: XNORM = MAX( ABS( XR1 )+ABS( XI1 ), ABS( XR2 )+ABS( XI2 ) )
487: *
488: * Further scaling if norm(A) norm(X) > overflow
489: *
490: IF( XNORM.GT.ONE .AND. CMAX.GT.ONE ) THEN
491: IF( XNORM.GT.BIGNUM / CMAX ) THEN
492: TEMP = CMAX / BIGNUM
493: X( 1, 1 ) = TEMP*X( 1, 1 )
494: X( 2, 1 ) = TEMP*X( 2, 1 )
495: X( 1, 2 ) = TEMP*X( 1, 2 )
496: X( 2, 2 ) = TEMP*X( 2, 2 )
497: XNORM = TEMP*XNORM
498: SCALE = TEMP*SCALE
499: END IF
500: END IF
501: END IF
502: END IF
503: *
504: RETURN
505: *
506: * End of DLALN2
507: *
508: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaln2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack