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: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       LOGICAL            LTRANS
011:       INTEGER            INFO, LDA, LDB, LDX, NA, NW
012:       DOUBLE PRECISION   CA, D1, D2, SCALE, SMIN, WI, WR, XNORM
013: *     ..
014: *     .. Array Arguments ..
015:       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), X( LDX, * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  DLALN2 solves a system of the form  (ca A - w D ) X = s B
022: *  or (ca A' - w D) X = s B   with possible scaling ("s") and
023: *  perturbation of A.  (A' means A-transpose.)
024: *
025: *  A is an NA x NA real matrix, ca is a real scalar, D is an NA x NA
026: *  real diagonal matrix, w is a real or complex value, and X and B are
027: *  NA x 1 matrices -- real if w is real, complex if w is complex.  NA
028: *  may be 1 or 2.
029: *
030: *  If w is complex, X and B are represented as NA x 2 matrices,
031: *  the first column of each being the real part and the second
032: *  being the imaginary part.
033: *
034: *  "s" is a scaling factor (.LE. 1), computed by DLALN2, which is
035: *  so chosen that X can be computed without overflow.  X is further
036: *  scaled if necessary to assure that norm(ca A - w D)*norm(X) is less
037: *  than overflow.
038: *
039: *  If both singular values of (ca A - w D) are less than SMIN,
040: *  SMIN*identity will be used instead of (ca A - w D).  If only one
041: *  singular value is less than SMIN, one element of (ca A - w D) will be
042: *  perturbed enough to make the smallest singular value roughly SMIN.
043: *  If both singular values are at least SMIN, (ca A - w D) will not be
044: *  perturbed.  In any case, the perturbation will be at most some small
045: *  multiple of max( SMIN, ulp*norm(ca A - w D) ).  The singular values
046: *  are computed by infinity-norm approximations, and thus will only be
047: *  correct to a factor of 2 or so.
048: *
049: *  Note: all input quantities are assumed to be smaller than overflow
050: *  by a reasonable factor.  (See BIGNUM.)
051: *
052: *  Arguments
053: *  ==========
054: *
055: *  LTRANS  (input) LOGICAL
056: *          =.TRUE.:  A-transpose will be used.
057: *          =.FALSE.: A will be used (not transposed.)
058: *
059: *  NA      (input) INTEGER
060: *          The size of the matrix A.  It may (only) be 1 or 2.
061: *
062: *  NW      (input) INTEGER
063: *          1 if "w" is real, 2 if "w" is complex.  It may only be 1
064: *          or 2.
065: *
066: *  SMIN    (input) DOUBLE PRECISION
067: *          The desired lower bound on the singular values of A.  This
068: *          should be a safe distance away from underflow or overflow,
069: *          say, between (underflow/machine precision) and  (machine
070: *          precision * overflow ).  (See BIGNUM and ULP.)
071: *
072: *  CA      (input) DOUBLE PRECISION
073: *          The coefficient c, which A is multiplied by.
074: *
075: *  A       (input) DOUBLE PRECISION array, dimension (LDA,NA)
076: *          The NA x NA matrix A.
077: *
078: *  LDA     (input) INTEGER
079: *          The leading dimension of A.  It must be at least NA.
080: *
081: *  D1      (input) DOUBLE PRECISION
082: *          The 1,1 element in the diagonal matrix D.
083: *
084: *  D2      (input) DOUBLE PRECISION
085: *          The 2,2 element in the diagonal matrix D.  Not used if NW=1.
086: *
087: *  B       (input) DOUBLE PRECISION array, dimension (LDB,NW)
088: *          The NA x NW matrix B (right-hand side).  If NW=2 ("w" is
089: *          complex), column 1 contains the real part of B and column 2
090: *          contains the imaginary part.
091: *
092: *  LDB     (input) INTEGER
093: *          The leading dimension of B.  It must be at least NA.
094: *
095: *  WR      (input) DOUBLE PRECISION
096: *          The real part of the scalar "w".
097: *
098: *  WI      (input) DOUBLE PRECISION
099: *          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
509: