LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ dlaln2()

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 (<= 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 NA=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.

Definition at line 216 of file dlaln2.f.

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