 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.

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.```

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 ..
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.