001:       SUBROUTINE DLALN2( LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B,
002:      $                   LDB, WR, WI, X, LDX, SCALE, XNORM, INFO )
003: *
004: *  -- LAPACK auxiliary routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       LOGICAL            LTRANS
010:       INTEGER            INFO, LDA, LDB, LDX, NA, NW
011:       DOUBLE PRECISION   CA, D1, D2, SCALE, SMIN, WI, WR, XNORM
012: *     ..
013: *     .. Array Arguments ..
014:       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), X( LDX, * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  DLALN2 solves a system of the form  (ca A - w D ) X = s B
021: *  or (ca A' - w D) X = s B   with possible scaling ("s") and
022: *  perturbation of A.  (A' means A-transpose.)
023: *
024: *  A is an NA x NA real matrix, ca is a real scalar, D is an NA x NA
025: *  real diagonal matrix, w is a real or complex value, and X and B are
026: *  NA x 1 matrices -- real if w is real, complex if w is complex.  NA
027: *  may be 1 or 2.
028: *
029: *  If w is complex, X and B are represented as NA x 2 matrices,
030: *  the first column of each being the real part and the second
031: *  being the imaginary part.
032: *
033: *  "s" is a scaling factor (.LE. 1), computed by DLALN2, which is
034: *  so chosen that X can be computed without overflow.  X is further
035: *  scaled if necessary to assure that norm(ca A - w D)*norm(X) is less
036: *  than overflow.
037: *
038: *  If both singular values of (ca A - w D) are less than SMIN,
039: *  SMIN*identity will be used instead of (ca A - w D).  If only one
040: *  singular value is less than SMIN, one element of (ca A - w D) will be
041: *  perturbed enough to make the smallest singular value roughly SMIN.
042: *  If both singular values are at least SMIN, (ca A - w D) will not be
043: *  perturbed.  In any case, the perturbation will be at most some small
044: *  multiple of max( SMIN, ulp*norm(ca A - w D) ).  The singular values
045: *  are computed by infinity-norm approximations, and thus will only be
046: *  correct to a factor of 2 or so.
047: *
048: *  Note: all input quantities are assumed to be smaller than overflow
049: *  by a reasonable factor.  (See BIGNUM.)
050: *
051: *  Arguments
052: *  ==========
053: *
054: *  LTRANS  (input) LOGICAL
055: *          =.TRUE.:  A-transpose will be used.
056: *          =.FALSE.: A will be used (not transposed.)
057: *
058: *  NA      (input) INTEGER
059: *          The size of the matrix A.  It may (only) be 1 or 2.
060: *
061: *  NW      (input) INTEGER
062: *          1 if "w" is real, 2 if "w" is complex.  It may only be 1
063: *          or 2.
064: *
065: *  SMIN    (input) DOUBLE PRECISION
066: *          The desired lower bound on the singular values of A.  This
067: *          should be a safe distance away from underflow or overflow,
068: *          say, between (underflow/machine precision) and  (machine
069: *          precision * overflow ).  (See BIGNUM and ULP.)
070: *
071: *  CA      (input) DOUBLE PRECISION
072: *          The coefficient c, which A is multiplied by.
073: *
074: *  A       (input) DOUBLE PRECISION array, dimension (LDA,NA)
075: *          The NA x NA matrix A.
076: *
077: *  LDA     (input) INTEGER
078: *          The leading dimension of A.  It must be at least NA.
079: *
080: *  D1      (input) DOUBLE PRECISION
081: *          The 1,1 element in the diagonal matrix D.
082: *
083: *  D2      (input) DOUBLE PRECISION
084: *          The 2,2 element in the diagonal matrix D.  Not used if NW=1.
085: *
086: *  B       (input) DOUBLE PRECISION array, dimension (LDB,NW)
087: *          The NA x NW matrix B (right-hand side).  If NW=2 ("w" is
088: *          complex), column 1 contains the real part of B and column 2
089: *          contains the imaginary part.
090: *
091: *  LDB     (input) INTEGER
092: *          The leading dimension of B.  It must be at least NA.
093: *
094: *  WR      (input) DOUBLE PRECISION
095: *          The real part of the scalar "w".
096: *
097: *  WI      (input) DOUBLE PRECISION
098: *          The imaginary part of the scalar "w".  Not used if NW=1.
099: *
100: *  X       (output) DOUBLE PRECISION array, dimension (LDX,NW)
101: *          The NA x NW matrix X (unknowns), as computed by DLALN2.
102: *          If NW=2 ("w" is complex), on exit, column 1 will contain
103: *          the real part of X and column 2 will contain the imaginary
104: *          part.
105: *
106: *  LDX     (input) INTEGER
107: *          The leading dimension of X.  It must be at least NA.
108: *
109: *  SCALE   (output) DOUBLE PRECISION
110: *          The scale factor that B must be multiplied by to insure
111: *          that overflow does not occur when computing X.  Thus,
112: *          (ca A - w D) X  will be SCALE*B, not B (ignoring
113: *          perturbations of A.)  It will be at most 1.
114: *
115: *  XNORM   (output) DOUBLE PRECISION
116: *          The infinity-norm of X, when X is regarded as an NA x NW
117: *          real matrix.
118: *
119: *  INFO    (output) INTEGER
120: *          An error flag.  It will be set to zero if no error occurs,
121: *          a negative number if an argument is in error, or a positive
122: *          number if  ca A - w D  had to be perturbed.
123: *          The possible values are:
124: *          = 0: No error occurred, and (ca A - w D) did not have to be
125: *                 perturbed.
126: *          = 1: (ca A - w D) had to be perturbed to make its smallest
127: *               (or only) singular value greater than SMIN.
128: *          NOTE: In the interests of speed, this routine does not
129: *                check the inputs for errors.
130: *
131: * =====================================================================
132: *
133: *     .. Parameters ..
134:       DOUBLE PRECISION   ZERO, ONE
135:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
136:       DOUBLE PRECISION   TWO
137:       PARAMETER          ( TWO = 2.0D0 )
138: *     ..
139: *     .. Local Scalars ..
140:       INTEGER            ICMAX, J
141:       DOUBLE PRECISION   BBND, BI1, BI2, BIGNUM, BNORM, BR1, BR2, CI21,
142:      $                   CI22, CMAX, CNORM, CR21, CR22, CSI, CSR, LI21,
143:      $                   LR21, SMINI, SMLNUM, TEMP, U22ABS, UI11, UI11R,
144:      $                   UI12, UI12S, UI22, UR11, UR11R, UR12, UR12S,
145:      $                   UR22, XI1, XI2, XR1, XR2
146: *     ..
147: *     .. Local Arrays ..
148:       LOGICAL            RSWAP( 4 ), ZSWAP( 4 )
149:       INTEGER            IPIVOT( 4, 4 )
150:       DOUBLE PRECISION   CI( 2, 2 ), CIV( 4 ), CR( 2, 2 ), CRV( 4 )
151: *     ..
152: *     .. External Functions ..
153:       DOUBLE PRECISION   DLAMCH
154:       EXTERNAL           DLAMCH
155: *     ..
156: *     .. External Subroutines ..
157:       EXTERNAL           DLADIV
158: *     ..
159: *     .. Intrinsic Functions ..
160:       INTRINSIC          ABS, MAX
161: *     ..
162: *     .. Equivalences ..
163:       EQUIVALENCE        ( CI( 1, 1 ), CIV( 1 ) ),
164:      $                   ( CR( 1, 1 ), CRV( 1 ) )
165: *     ..
166: *     .. Data statements ..
167:       DATA               ZSWAP / .FALSE., .FALSE., .TRUE., .TRUE. /
168:       DATA               RSWAP / .FALSE., .TRUE., .FALSE., .TRUE. /
169:       DATA               IPIVOT / 1, 2, 3, 4, 2, 1, 4, 3, 3, 4, 1, 2, 4,
170:      $                   3, 2, 1 /
171: *     ..
172: *     .. Executable Statements ..
173: *
174: *     Compute BIGNUM
175: *
176:       SMLNUM = TWO*DLAMCH( 'Safe minimum' )
177:       BIGNUM = ONE / SMLNUM
178:       SMINI = MAX( SMIN, SMLNUM )
179: *
180: *     Don't check for input errors
181: *
182:       INFO = 0
183: *
184: *     Standard Initializations
185: *
186:       SCALE = ONE
187: *
188:       IF( NA.EQ.1 ) THEN
189: *
190: *        1 x 1  (i.e., scalar) system   C X = B
191: *
192:          IF( NW.EQ.1 ) THEN
193: *
194: *           Real 1x1 system.
195: *
196: *           C = ca A - w D
197: *
198:             CSR = CA*A( 1, 1 ) - WR*D1
199:             CNORM = ABS( CSR )
200: *
201: *           If | C | < SMINI, use C = SMINI
202: *
203:             IF( CNORM.LT.SMINI ) THEN
204:                CSR = SMINI
205:                CNORM = SMINI
206:                INFO = 1
207:             END IF
208: *
209: *           Check scaling for  X = B / C
210: *
211:             BNORM = ABS( B( 1, 1 ) )
212:             IF( CNORM.LT.ONE .AND. BNORM.GT.ONE ) THEN
213:                IF( BNORM.GT.BIGNUM*CNORM )
214:      $            SCALE = ONE / BNORM
215:             END IF
216: *
217: *           Compute X
218: *
219:             X( 1, 1 ) = ( B( 1, 1 )*SCALE ) / CSR
220:             XNORM = ABS( X( 1, 1 ) )
221:          ELSE
222: *
223: *           Complex 1x1 system (w is complex)
224: *
225: *           C = ca A - w D
226: *
227:             CSR = CA*A( 1, 1 ) - WR*D1
228:             CSI = -WI*D1
229:             CNORM = ABS( CSR ) + ABS( CSI )
230: *
231: *           If | C | < SMINI, use C = SMINI
232: *
233:             IF( CNORM.LT.SMINI ) THEN
234:                CSR = SMINI
235:                CSI = ZERO
236:                CNORM = SMINI
237:                INFO = 1
238:             END IF
239: *
240: *           Check scaling for  X = B / C
241: *
242:             BNORM = ABS( B( 1, 1 ) ) + ABS( B( 1, 2 ) )
243:             IF( CNORM.LT.ONE .AND. BNORM.GT.ONE ) THEN
244:                IF( BNORM.GT.BIGNUM*CNORM )
245:      $            SCALE = ONE / BNORM
246:             END IF
247: *
248: *           Compute X
249: *
250:             CALL DLADIV( SCALE*B( 1, 1 ), SCALE*B( 1, 2 ), CSR, CSI,
251:      $                   X( 1, 1 ), X( 1, 2 ) )
252:             XNORM = ABS( X( 1, 1 ) ) + ABS( X( 1, 2 ) )
253:          END IF
254: *
255:       ELSE
256: *
257: *        2x2 System
258: *
259: *        Compute the real part of  C = ca A - w D  (or  ca A' - w D )
260: *
261:          CR( 1, 1 ) = CA*A( 1, 1 ) - WR*D1
262:          CR( 2, 2 ) = CA*A( 2, 2 ) - WR*D2
263:          IF( LTRANS ) THEN
264:             CR( 1, 2 ) = CA*A( 2, 1 )
265:             CR( 2, 1 ) = CA*A( 1, 2 )
266:          ELSE
267:             CR( 2, 1 ) = CA*A( 2, 1 )
268:             CR( 1, 2 ) = CA*A( 1, 2 )
269:          END IF
270: *
271:          IF( NW.EQ.1 ) THEN
272: *
273: *           Real 2x2 system  (w is real)
274: *
275: *           Find the largest element in C
276: *
277:             CMAX = ZERO
278:             ICMAX = 0
279: *
280:             DO 10 J = 1, 4
281:                IF( ABS( CRV( J ) ).GT.CMAX ) THEN
282:                   CMAX = ABS( CRV( J ) )
283:                   ICMAX = J
284:                END IF
285:    10       CONTINUE
286: *
287: *           If norm(C) < SMINI, use SMINI*identity.
288: *
289:             IF( CMAX.LT.SMINI ) THEN
290:                BNORM = MAX( ABS( B( 1, 1 ) ), ABS( B( 2, 1 ) ) )
291:                IF( SMINI.LT.ONE .AND. BNORM.GT.ONE ) THEN
292:                   IF( BNORM.GT.BIGNUM*SMINI )
293:      $               SCALE = ONE / BNORM
294:                END IF
295:                TEMP = SCALE / SMINI
296:                X( 1, 1 ) = TEMP*B( 1, 1 )
297:                X( 2, 1 ) = TEMP*B( 2, 1 )
298:                XNORM = TEMP*BNORM
299:                INFO = 1
300:                RETURN
301:             END IF
302: *
303: *           Gaussian elimination with complete pivoting.
304: *
305:             UR11 = CRV( ICMAX )
306:             CR21 = CRV( IPIVOT( 2, ICMAX ) )
307:             UR12 = CRV( IPIVOT( 3, ICMAX ) )
308:             CR22 = CRV( IPIVOT( 4, ICMAX ) )
309:             UR11R = ONE / UR11
310:             LR21 = UR11R*CR21
311:             UR22 = CR22 - UR12*LR21
312: *
313: *           If smaller pivot < SMINI, use SMINI
314: *
315:             IF( ABS( UR22 ).LT.SMINI ) THEN
316:                UR22 = SMINI
317:                INFO = 1
318:             END IF
319:             IF( RSWAP( ICMAX ) ) THEN
320:                BR1 = B( 2, 1 )
321:                BR2 = B( 1, 1 )
322:             ELSE
323:                BR1 = B( 1, 1 )
324:                BR2 = B( 2, 1 )
325:             END IF
326:             BR2 = BR2 - LR21*BR1
327:             BBND = MAX( ABS( BR1*( UR22*UR11R ) ), ABS( BR2 ) )
328:             IF( BBND.GT.ONE .AND. ABS( UR22 ).LT.ONE ) THEN
329:                IF( BBND.GE.BIGNUM*ABS( UR22 ) )
330:      $            SCALE = ONE / BBND
331:             END IF
332: *
333:             XR2 = ( BR2*SCALE ) / UR22
334:             XR1 = ( SCALE*BR1 )*UR11R - XR2*( UR11R*UR12 )
335:             IF( ZSWAP( ICMAX ) ) THEN
336:                X( 1, 1 ) = XR2
337:                X( 2, 1 ) = XR1
338:             ELSE
339:                X( 1, 1 ) = XR1
340:                X( 2, 1 ) = XR2
341:             END IF
342:             XNORM = MAX( ABS( XR1 ), ABS( XR2 ) )
343: *
344: *           Further scaling if  norm(A) norm(X) > overflow
345: *
346:             IF( XNORM.GT.ONE .AND. CMAX.GT.ONE ) THEN
347:                IF( XNORM.GT.BIGNUM / CMAX ) THEN
348:                   TEMP = CMAX / BIGNUM
349:                   X( 1, 1 ) = TEMP*X( 1, 1 )
350:                   X( 2, 1 ) = TEMP*X( 2, 1 )
351:                   XNORM = TEMP*XNORM
352:                   SCALE = TEMP*SCALE
353:                END IF
354:             END IF
355:          ELSE
356: *
357: *           Complex 2x2 system  (w is complex)
358: *
359: *           Find the largest element in C
360: *
361:             CI( 1, 1 ) = -WI*D1
362:             CI( 2, 1 ) = ZERO
363:             CI( 1, 2 ) = ZERO
364:             CI( 2, 2 ) = -WI*D2
365:             CMAX = ZERO
366:             ICMAX = 0
367: *
368:             DO 20 J = 1, 4
369:                IF( ABS( CRV( J ) )+ABS( CIV( J ) ).GT.CMAX ) THEN
370:                   CMAX = ABS( CRV( J ) ) + ABS( CIV( J ) )
371:                   ICMAX = J
372:                END IF
373:    20       CONTINUE
374: *
375: *           If norm(C) < SMINI, use SMINI*identity.
376: *
377:             IF( CMAX.LT.SMINI ) THEN
378:                BNORM = MAX( ABS( B( 1, 1 ) )+ABS( B( 1, 2 ) ),
379:      $                 ABS( B( 2, 1 ) )+ABS( B( 2, 2 ) ) )
380:                IF( SMINI.LT.ONE .AND. BNORM.GT.ONE ) THEN
381:                   IF( BNORM.GT.BIGNUM*SMINI )
382:      $               SCALE = ONE / BNORM
383:                END IF
384:                TEMP = SCALE / SMINI
385:                X( 1, 1 ) = TEMP*B( 1, 1 )
386:                X( 2, 1 ) = TEMP*B( 2, 1 )
387:                X( 1, 2 ) = TEMP*B( 1, 2 )
388:                X( 2, 2 ) = TEMP*B( 2, 2 )
389:                XNORM = TEMP*BNORM
390:                INFO = 1
391:                RETURN
392:             END IF
393: *
394: *           Gaussian elimination with complete pivoting.
395: *
396:             UR11 = CRV( ICMAX )
397:             UI11 = CIV( ICMAX )
398:             CR21 = CRV( IPIVOT( 2, ICMAX ) )
399:             CI21 = CIV( IPIVOT( 2, ICMAX ) )
400:             UR12 = CRV( IPIVOT( 3, ICMAX ) )
401:             UI12 = CIV( IPIVOT( 3, ICMAX ) )
402:             CR22 = CRV( IPIVOT( 4, ICMAX ) )
403:             CI22 = CIV( IPIVOT( 4, ICMAX ) )
404:             IF( ICMAX.EQ.1 .OR. ICMAX.EQ.4 ) THEN
405: *
406: *              Code when off-diagonals of pivoted C are real
407: *
408:                IF( ABS( UR11 ).GT.ABS( UI11 ) ) THEN
409:                   TEMP = UI11 / UR11
410:                   UR11R = ONE / ( UR11*( ONE+TEMP**2 ) )
411:                   UI11R = -TEMP*UR11R
412:                ELSE
413:                   TEMP = UR11 / UI11
414:                   UI11R = -ONE / ( UI11*( ONE+TEMP**2 ) )
415:                   UR11R = -TEMP*UI11R
416:                END IF
417:                LR21 = CR21*UR11R
418:                LI21 = CR21*UI11R
419:                UR12S = UR12*UR11R
420:                UI12S = UR12*UI11R
421:                UR22 = CR22 - UR12*LR21
422:                UI22 = CI22 - UR12*LI21
423:             ELSE
424: *
425: *              Code when diagonals of pivoted C are real
426: *
427:                UR11R = ONE / UR11
428:                UI11R = ZERO
429:                LR21 = CR21*UR11R
430:                LI21 = CI21*UR11R
431:                UR12S = UR12*UR11R
432:                UI12S = UI12*UR11R
433:                UR22 = CR22 - UR12*LR21 + UI12*LI21
434:                UI22 = -UR12*LI21 - UI12*LR21
435:             END IF
436:             U22ABS = ABS( UR22 ) + ABS( UI22 )
437: *
438: *           If smaller pivot < SMINI, use SMINI
439: *
440:             IF( U22ABS.LT.SMINI ) THEN
441:                UR22 = SMINI
442:                UI22 = ZERO
443:                INFO = 1
444:             END IF
445:             IF( RSWAP( ICMAX ) ) THEN
446:                BR2 = B( 1, 1 )
447:                BR1 = B( 2, 1 )
448:                BI2 = B( 1, 2 )
449:                BI1 = B( 2, 2 )
450:             ELSE
451:                BR1 = B( 1, 1 )
452:                BR2 = B( 2, 1 )
453:                BI1 = B( 1, 2 )
454:                BI2 = B( 2, 2 )
455:             END IF
456:             BR2 = BR2 - LR21*BR1 + LI21*BI1
457:             BI2 = BI2 - LI21*BR1 - LR21*BI1
458:             BBND = MAX( ( ABS( BR1 )+ABS( BI1 ) )*
459:      $             ( U22ABS*( ABS( UR11R )+ABS( UI11R ) ) ),
460:      $             ABS( BR2 )+ABS( BI2 ) )
461:             IF( BBND.GT.ONE .AND. U22ABS.LT.ONE ) THEN
462:                IF( BBND.GE.BIGNUM*U22ABS ) THEN
463:                   SCALE = ONE / BBND
464:                   BR1 = SCALE*BR1
465:                   BI1 = SCALE*BI1
466:                   BR2 = SCALE*BR2
467:                   BI2 = SCALE*BI2
468:                END IF
469:             END IF
470: *
471:             CALL DLADIV( BR2, BI2, UR22, UI22, XR2, XI2 )
472:             XR1 = UR11R*BR1 - UI11R*BI1 - UR12S*XR2 + UI12S*XI2
473:             XI1 = UI11R*BR1 + UR11R*BI1 - UI12S*XR2 - UR12S*XI2
474:             IF( ZSWAP( ICMAX ) ) THEN
475:                X( 1, 1 ) = XR2
476:                X( 2, 1 ) = XR1
477:                X( 1, 2 ) = XI2
478:                X( 2, 2 ) = XI1
479:             ELSE
480:                X( 1, 1 ) = XR1
481:                X( 2, 1 ) = XR2
482:                X( 1, 2 ) = XI1
483:                X( 2, 2 ) = XI2
484:             END IF
485:             XNORM = MAX( ABS( XR1 )+ABS( XI1 ), ABS( XR2 )+ABS( XI2 ) )
486: *
487: *           Further scaling if  norm(A) norm(X) > overflow
488: *
489:             IF( XNORM.GT.ONE .AND. CMAX.GT.ONE ) THEN
490:                IF( XNORM.GT.BIGNUM / CMAX ) THEN
491:                   TEMP = CMAX / BIGNUM
492:                   X( 1, 1 ) = TEMP*X( 1, 1 )
493:                   X( 2, 1 ) = TEMP*X( 2, 1 )
494:                   X( 1, 2 ) = TEMP*X( 1, 2 )
495:                   X( 2, 2 ) = TEMP*X( 2, 2 )
496:                   XNORM = TEMP*XNORM
497:                   SCALE = TEMP*SCALE
498:                END IF
499:             END IF
500:          END IF
501:       END IF
502: *
503:       RETURN
504: *
505: *     End of DLALN2
506: *
507:       END
508: