LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dlaln2 ( logical  LTRANS,
integer  NA,
integer  NW,
double precision  SMIN,
double precision  CA,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision  D1,
double precision  D2,
double precision, dimension( ldb, * )  B,
integer  LDB,
double precision  WR,
double precision  WI,
double precision, dimension( ldx, * )  X,
integer  LDX,
double precision  SCALE,
double precision  XNORM,
integer  INFO 
)

DLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form.

Download DLALN2 + dependencies [TGZ] [ZIP] [TXT]

Purpose: 
 DLALN2 solves a system of the form  (ca A - w D ) X = s B
 or (ca A**T - w D) X = s B   with possible scaling ("s") and
 perturbation of A.  (A**T means A-transpose.)

 A is an NA x NA real matrix, ca is a real scalar, D is an NA x NA
 real diagonal matrix, w is a real or complex value, and X and B are
 NA x 1 matrices -- real if w is real, complex if w is complex.  NA
 may be 1 or 2.

 If w is complex, X and B are represented as NA x 2 matrices,
 the first column of each being the real part and the second
 being the imaginary part.

 "s" is a scaling factor (.LE. 1), computed by DLALN2, which is
 so chosen that X can be computed without overflow.  X is further
 scaled if necessary to assure that norm(ca A - w D)*norm(X) is less
 than overflow.

 If both singular values of (ca A - w D) are less than SMIN,
 SMIN*identity will be used instead of (ca A - w D).  If only one
 singular value is less than SMIN, one element of (ca A - w D) will be
 perturbed enough to make the smallest singular value roughly SMIN.
 If both singular values are at least SMIN, (ca A - w D) will not be
 perturbed.  In any case, the perturbation will be at most some small
 multiple of max( SMIN, ulp*norm(ca A - w D) ).  The singular values
 are computed by infinity-norm approximations, and thus will only be
 correct to a factor of 2 or so.

 Note: all input quantities are assumed to be smaller than overflow
 by a reasonable factor.  (See BIGNUM.)
Parameters
[in]LTRANS
          LTRANS is LOGICAL
          =.TRUE.:  A-transpose will be used.
          =.FALSE.: A will be used (not transposed.)
[in]NA
          NA is INTEGER
          The size of the matrix A.  It may (only) be 1 or 2.
[in]NW
          NW is INTEGER
          1 if "w" is real, 2 if "w" is complex.  It may only be 1
          or 2.
[in]SMIN
          SMIN is DOUBLE PRECISION
          The desired lower bound on the singular values of A.  This
          should be a safe distance away from underflow or overflow,
          say, between (underflow/machine precision) and  (machine
          precision * overflow ).  (See BIGNUM and ULP.)
[in]CA
          CA is DOUBLE PRECISION
          The coefficient c, which A is multiplied by.
[in]A
          A is DOUBLE PRECISION array, dimension (LDA,NA)
          The NA x NA matrix A.
[in]LDA
          LDA is INTEGER
          The leading dimension of A.  It must be at least NA.
[in]D1
          D1 is DOUBLE PRECISION
          The 1,1 element in the diagonal matrix D.
[in]D2
          D2 is DOUBLE PRECISION
          The 2,2 element in the diagonal matrix D.  Not used if NW=1.
[in]B
          B is DOUBLE PRECISION array, dimension (LDB,NW)
          The NA x NW matrix B (right-hand side).  If NW=2 ("w" is
          complex), column 1 contains the real part of B and column 2
          contains the imaginary part.
[in]LDB
          LDB is INTEGER
          The leading dimension of B.  It must be at least NA.
[in]WR
          WR is DOUBLE PRECISION
          The real part of the scalar "w".
[in]WI
          WI is DOUBLE PRECISION
          The imaginary part of the scalar "w".  Not used if NW=1.
[out]X
          X is DOUBLE PRECISION array, dimension (LDX,NW)
          The NA x NW matrix X (unknowns), as computed by DLALN2.
          If NW=2 ("w" is complex), on exit, column 1 will contain
          the real part of X and column 2 will contain the imaginary
          part.
[in]LDX
          LDX is INTEGER
          The leading dimension of X.  It must be at least NA.
[out]SCALE
          SCALE is DOUBLE PRECISION
          The scale factor that B must be multiplied by to insure
          that overflow does not occur when computing X.  Thus,
          (ca A - w D) X  will be SCALE*B, not B (ignoring
          perturbations of A.)  It will be at most 1.
[out]XNORM
          XNORM is DOUBLE PRECISION
          The infinity-norm of X, when X is regarded as an NA x NW
          real matrix.
[out]INFO
          INFO is INTEGER
          An error flag.  It will be set to zero if no error occurs,
          a negative number if an argument is in error, or a positive
          number if  ca A - w D  had to be perturbed.
          The possible values are:
          = 0: No error occurred, and (ca A - w D) did not have to be
                 perturbed.
          = 1: (ca A - w D) had to be perturbed to make its smallest
               (or only) singular value greater than SMIN.
          NOTE: In the interests of speed, this routine does not
                check the inputs for errors.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 220 of file dlaln2.f.

220 *
221 * -- LAPACK auxiliary routine (version 3.4.2) --
222 * -- LAPACK is a software package provided by Univ. of Tennessee, --
223 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
224 * September 2012
225 *
226 * .. Scalar Arguments ..
227  LOGICAL ltrans
228  INTEGER info, lda, ldb, ldx, na, nw
229  DOUBLE PRECISION ca, d1, d2, scale, smin, wi, wr, xnorm
230 * ..
231 * .. Array Arguments ..
232  DOUBLE PRECISION a( lda, * ), b( ldb, * ), x( ldx, * )
233 * ..
234 *
235 * =====================================================================
236 *
237 * .. Parameters ..
238  DOUBLE PRECISION zero, one
239  parameter ( zero = 0.0d0, one = 1.0d0 )
240  DOUBLE PRECISION two
241  parameter ( two = 2.0d0 )
242 * ..
243 * .. Local Scalars ..
244  INTEGER icmax, j
245  DOUBLE PRECISION bbnd, bi1, bi2, bignum, bnorm, br1, br2, ci21,
246  $ ci22, cmax, cnorm, cr21, cr22, csi, csr, li21,
247  $ lr21, smini, smlnum, temp, u22abs, ui11, ui11r,
248  $ ui12, ui12s, ui22, ur11, ur11r, ur12, ur12s,
249  $ ur22, xi1, xi2, xr1, xr2
250 * ..
251 * .. Local Arrays ..
252  LOGICAL rswap( 4 ), zswap( 4 )
253  INTEGER ipivot( 4, 4 )
254  DOUBLE PRECISION ci( 2, 2 ), civ( 4 ), cr( 2, 2 ), crv( 4 )
255 * ..
256 * .. External Functions ..
257  DOUBLE PRECISION dlamch
258  EXTERNAL dlamch
259 * ..
260 * .. External Subroutines ..
261  EXTERNAL dladiv
262 * ..
263 * .. Intrinsic Functions ..
264  INTRINSIC abs, max
265 * ..
266 * .. Equivalences ..
267  equivalence ( ci( 1, 1 ), civ( 1 ) ),
268  $ ( cr( 1, 1 ), crv( 1 ) )
269 * ..
270 * .. Data statements ..
271  DATA zswap / .false., .false., .true., .true. /
272  DATA rswap / .false., .true., .false., .true. /
273  DATA ipivot / 1, 2, 3, 4, 2, 1, 4, 3, 3, 4, 1, 2, 4,
274  $ 3, 2, 1 /
275 * ..
276 * .. Executable Statements ..
277 *
278 * Compute BIGNUM
279 *
280  smlnum = two*dlamch( 'Safe minimum' )
281  bignum = one / smlnum
282  smini = max( smin, smlnum )
283 *
284 * Don't check for input errors
285 *
286  info = 0
287 *
288 * Standard Initializations
289 *
290  scale = one
291 *
292  IF( na.EQ.1 ) THEN
293 *
294 * 1 x 1 (i.e., scalar) system C X = B
295 *
296  IF( nw.EQ.1 ) THEN
297 *
298 * Real 1x1 system.
299 *
300 * C = ca A - w D
301 *
302  csr = ca*a( 1, 1 ) - wr*d1
303  cnorm = abs( csr )
304 *
305 * If | C | < SMINI, use C = SMINI
306 *
307  IF( cnorm.LT.smini ) THEN
308  csr = smini
309  cnorm = smini
310  info = 1
311  END IF
312 *
313 * Check scaling for X = B / C
314 *
315  bnorm = abs( b( 1, 1 ) )
316  IF( cnorm.LT.one .AND. bnorm.GT.one ) THEN
317  IF( bnorm.GT.bignum*cnorm )
318  $ scale = one / bnorm
319  END IF
320 *
321 * Compute X
322 *
323  x( 1, 1 ) = ( b( 1, 1 )*scale ) / csr
324  xnorm = abs( x( 1, 1 ) )
325  ELSE
326 *
327 * Complex 1x1 system (w is complex)
328 *
329 * C = ca A - w D
330 *
331  csr = ca*a( 1, 1 ) - wr*d1
332  csi = -wi*d1
333  cnorm = abs( csr ) + abs( csi )
334 *
335 * If | C | < SMINI, use C = SMINI
336 *
337  IF( cnorm.LT.smini ) THEN
338  csr = smini
339  csi = zero
340  cnorm = smini
341  info = 1
342  END IF
343 *
344 * Check scaling for X = B / C
345 *
346  bnorm = abs( b( 1, 1 ) ) + abs( b( 1, 2 ) )
347  IF( cnorm.LT.one .AND. bnorm.GT.one ) THEN
348  IF( bnorm.GT.bignum*cnorm )
349  $ scale = one / bnorm
350  END IF
351 *
352 * Compute X
353 *
354  CALL dladiv( scale*b( 1, 1 ), scale*b( 1, 2 ), csr, csi,
355  $ x( 1, 1 ), x( 1, 2 ) )
356  xnorm = abs( x( 1, 1 ) ) + abs( x( 1, 2 ) )
357  END IF
358 *
359  ELSE
360 *
361 * 2x2 System
362 *
363 * Compute the real part of C = ca A - w D (or ca A**T - w D )
364 *
365  cr( 1, 1 ) = ca*a( 1, 1 ) - wr*d1
366  cr( 2, 2 ) = ca*a( 2, 2 ) - wr*d2
367  IF( ltrans ) THEN
368  cr( 1, 2 ) = ca*a( 2, 1 )
369  cr( 2, 1 ) = ca*a( 1, 2 )
370  ELSE
371  cr( 2, 1 ) = ca*a( 2, 1 )
372  cr( 1, 2 ) = ca*a( 1, 2 )
373  END IF
374 *
375  IF( nw.EQ.1 ) THEN
376 *
377 * Real 2x2 system (w is real)
378 *
379 * Find the largest element in C
380 *
381  cmax = zero
382  icmax = 0
383 *
384  DO 10 j = 1, 4
385  IF( abs( crv( j ) ).GT.cmax ) THEN
386  cmax = abs( crv( j ) )
387  icmax = j
388  END IF
389  10 CONTINUE
390 *
391 * If norm(C) < SMINI, use SMINI*identity.
392 *
393  IF( cmax.LT.smini ) THEN
394  bnorm = max( abs( b( 1, 1 ) ), abs( b( 2, 1 ) ) )
395  IF( smini.LT.one .AND. bnorm.GT.one ) THEN
396  IF( bnorm.GT.bignum*smini )
397  $ scale = one / bnorm
398  END IF
399  temp = scale / smini
400  x( 1, 1 ) = temp*b( 1, 1 )
401  x( 2, 1 ) = temp*b( 2, 1 )
402  xnorm = temp*bnorm
403  info = 1
404  RETURN
405  END IF
406 *
407 * Gaussian elimination with complete pivoting.
408 *
409  ur11 = crv( icmax )
410  cr21 = crv( ipivot( 2, icmax ) )
411  ur12 = crv( ipivot( 3, icmax ) )
412  cr22 = crv( ipivot( 4, icmax ) )
413  ur11r = one / ur11
414  lr21 = ur11r*cr21
415  ur22 = cr22 - ur12*lr21
416 *
417 * If smaller pivot < SMINI, use SMINI
418 *
419  IF( abs( ur22 ).LT.smini ) THEN
420  ur22 = smini
421  info = 1
422  END IF
423  IF( rswap( icmax ) ) THEN
424  br1 = b( 2, 1 )
425  br2 = b( 1, 1 )
426  ELSE
427  br1 = b( 1, 1 )
428  br2 = b( 2, 1 )
429  END IF
430  br2 = br2 - lr21*br1
431  bbnd = max( abs( br1*( ur22*ur11r ) ), abs( br2 ) )
432  IF( bbnd.GT.one .AND. abs( ur22 ).LT.one ) THEN
433  IF( bbnd.GE.bignum*abs( ur22 ) )
434  $ scale = one / bbnd
435  END IF
436 *
437  xr2 = ( br2*scale ) / ur22
438  xr1 = ( scale*br1 )*ur11r - xr2*( ur11r*ur12 )
439  IF( zswap( icmax ) ) THEN
440  x( 1, 1 ) = xr2
441  x( 2, 1 ) = xr1
442  ELSE
443  x( 1, 1 ) = xr1
444  x( 2, 1 ) = xr2
445  END IF
446  xnorm = max( abs( xr1 ), abs( xr2 ) )
447 *
448 * Further scaling if norm(A) norm(X) > overflow
449 *
450  IF( xnorm.GT.one .AND. cmax.GT.one ) THEN
451  IF( xnorm.GT.bignum / cmax ) THEN
452  temp = cmax / bignum
453  x( 1, 1 ) = temp*x( 1, 1 )
454  x( 2, 1 ) = temp*x( 2, 1 )
455  xnorm = temp*xnorm
456  scale = temp*scale
457  END IF
458  END IF
459  ELSE
460 *
461 * Complex 2x2 system (w is complex)
462 *
463 * Find the largest element in C
464 *
465  ci( 1, 1 ) = -wi*d1
466  ci( 2, 1 ) = zero
467  ci( 1, 2 ) = zero
468  ci( 2, 2 ) = -wi*d2
469  cmax = zero
470  icmax = 0
471 *
472  DO 20 j = 1, 4
473  IF( abs( crv( j ) )+abs( civ( j ) ).GT.cmax ) THEN
474  cmax = abs( crv( j ) ) + abs( civ( j ) )
475  icmax = j
476  END IF
477  20 CONTINUE
478 *
479 * If norm(C) < SMINI, use SMINI*identity.
480 *
481  IF( cmax.LT.smini ) THEN
482  bnorm = max( abs( b( 1, 1 ) )+abs( b( 1, 2 ) ),
483  $ abs( b( 2, 1 ) )+abs( b( 2, 2 ) ) )
484  IF( smini.LT.one .AND. bnorm.GT.one ) THEN
485  IF( bnorm.GT.bignum*smini )
486  $ scale = one / bnorm
487  END IF
488  temp = scale / smini
489  x( 1, 1 ) = temp*b( 1, 1 )
490  x( 2, 1 ) = temp*b( 2, 1 )
491  x( 1, 2 ) = temp*b( 1, 2 )
492  x( 2, 2 ) = temp*b( 2, 2 )
493  xnorm = temp*bnorm
494  info = 1
495  RETURN
496  END IF
497 *
498 * Gaussian elimination with complete pivoting.
499 *
500  ur11 = crv( icmax )
501  ui11 = civ( icmax )
502  cr21 = crv( ipivot( 2, icmax ) )
503  ci21 = civ( ipivot( 2, icmax ) )
504  ur12 = crv( ipivot( 3, icmax ) )
505  ui12 = civ( ipivot( 3, icmax ) )
506  cr22 = crv( ipivot( 4, icmax ) )
507  ci22 = civ( ipivot( 4, icmax ) )
508  IF( icmax.EQ.1 .OR. icmax.EQ.4 ) THEN
509 *
510 * Code when off-diagonals of pivoted C are real
511 *
512  IF( abs( ur11 ).GT.abs( ui11 ) ) THEN
513  temp = ui11 / ur11
514  ur11r = one / ( ur11*( one+temp**2 ) )
515  ui11r = -temp*ur11r
516  ELSE
517  temp = ur11 / ui11
518  ui11r = -one / ( ui11*( one+temp**2 ) )
519  ur11r = -temp*ui11r
520  END IF
521  lr21 = cr21*ur11r
522  li21 = cr21*ui11r
523  ur12s = ur12*ur11r
524  ui12s = ur12*ui11r
525  ur22 = cr22 - ur12*lr21
526  ui22 = ci22 - ur12*li21
527  ELSE
528 *
529 * Code when diagonals of pivoted C are real
530 *
531  ur11r = one / ur11
532  ui11r = zero
533  lr21 = cr21*ur11r
534  li21 = ci21*ur11r
535  ur12s = ur12*ur11r
536  ui12s = ui12*ur11r
537  ur22 = cr22 - ur12*lr21 + ui12*li21
538  ui22 = -ur12*li21 - ui12*lr21
539  END IF
540  u22abs = abs( ur22 ) + abs( ui22 )
541 *
542 * If smaller pivot < SMINI, use SMINI
543 *
544  IF( u22abs.LT.smini ) THEN
545  ur22 = smini
546  ui22 = zero
547  info = 1
548  END IF
549  IF( rswap( icmax ) ) THEN
550  br2 = b( 1, 1 )
551  br1 = b( 2, 1 )
552  bi2 = b( 1, 2 )
553  bi1 = b( 2, 2 )
554  ELSE
555  br1 = b( 1, 1 )
556  br2 = b( 2, 1 )
557  bi1 = b( 1, 2 )
558  bi2 = b( 2, 2 )
559  END IF
560  br2 = br2 - lr21*br1 + li21*bi1
561  bi2 = bi2 - li21*br1 - lr21*bi1
562  bbnd = max( ( abs( br1 )+abs( bi1 ) )*
563  $ ( u22abs*( abs( ur11r )+abs( ui11r ) ) ),
564  $ abs( br2 )+abs( bi2 ) )
565  IF( bbnd.GT.one .AND. u22abs.LT.one ) THEN
566  IF( bbnd.GE.bignum*u22abs ) THEN
567  scale = one / bbnd
568  br1 = scale*br1
569  bi1 = scale*bi1
570  br2 = scale*br2
571  bi2 = scale*bi2
572  END IF
573  END IF
574 *
575  CALL dladiv( br2, bi2, ur22, ui22, xr2, xi2 )
576  xr1 = ur11r*br1 - ui11r*bi1 - ur12s*xr2 + ui12s*xi2
577  xi1 = ui11r*br1 + ur11r*bi1 - ui12s*xr2 - ur12s*xi2
578  IF( zswap( icmax ) ) THEN
579  x( 1, 1 ) = xr2
580  x( 2, 1 ) = xr1
581  x( 1, 2 ) = xi2
582  x( 2, 2 ) = xi1
583  ELSE
584  x( 1, 1 ) = xr1
585  x( 2, 1 ) = xr2
586  x( 1, 2 ) = xi1
587  x( 2, 2 ) = xi2
588  END IF
589  xnorm = max( abs( xr1 )+abs( xi1 ), abs( xr2 )+abs( xi2 ) )
590 *
591 * Further scaling if norm(A) norm(X) > overflow
592 *
593  IF( xnorm.GT.one .AND. cmax.GT.one ) THEN
594  IF( xnorm.GT.bignum / cmax ) THEN
595  temp = cmax / bignum
596  x( 1, 1 ) = temp*x( 1, 1 )
597  x( 2, 1 ) = temp*x( 2, 1 )
598  x( 1, 2 ) = temp*x( 1, 2 )
599  x( 2, 2 ) = temp*x( 2, 2 )
600  xnorm = temp*xnorm
601  scale = temp*scale
602  END IF
603  END IF
604  END IF
605  END IF
606 *
607  RETURN
608 *
609 * End of DLALN2
610 *
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
dladiv
subroutine dladiv(A, B, C, D, P, Q)
DLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
Definition: dladiv.f:93
zswap
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:52

Here is the call graph for this function:

Here is the caller graph for this function: