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

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

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

Purpose:
``` SLALN2 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 SLALN2, 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 REAL 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 REAL The coefficient c, which A is multiplied by.``` [in] A ``` A is REAL 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 REAL The 1,1 element in the diagonal matrix D.``` [in] D2 ``` D2 is REAL The 2,2 element in the diagonal matrix D. Not used if NW=1.``` [in] B ``` B is REAL 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 REAL The real part of the scalar "w".``` [in] WI ``` WI is REAL The imaginary part of the scalar "w". Not used if NW=1.``` [out] X ``` X is REAL array, dimension (LDX,NW) The NA x NW matrix X (unknowns), as computed by SLALN2. 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 REAL 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 REAL 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.```
Date
September 2012

Definition at line 220 of file slaln2.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  REAL ca, d1, d2, scale, smin, wi, wr, xnorm
230 * ..
231 * .. Array Arguments ..
232  REAL a( lda, * ), b( ldb, * ), x( ldx, * )
233 * ..
234 *
235 * =====================================================================
236 *
237 * .. Parameters ..
238  REAL zero, one
239  parameter ( zero = 0.0e0, one = 1.0e0 )
240  REAL two
241  parameter ( two = 2.0e0 )
242 * ..
243 * .. Local Scalars ..
244  INTEGER icmax, j
245  REAL 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 cswap( 4 ), rswap( 4 )
253  INTEGER ipivot( 4, 4 )
254  REAL ci( 2, 2 ), civ( 4 ), cr( 2, 2 ), crv( 4 )
255 * ..
256 * .. External Functions ..
257  REAL slamch
258  EXTERNAL slamch
259 * ..
260 * .. External Subroutines ..
261  EXTERNAL sladiv
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 cswap / .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*slamch( '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 sladiv( 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( cswap( 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 sladiv( 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( cswap( 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 SLALN2
610 *
subroutine sladiv(A, B, C, D, P, Q)
SLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
Definition: sladiv.f:93
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:52
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69

Here is the call graph for this function:

Here is the caller graph for this function: