KTH framework for Nek5000 toolboxes; testing version  0.0.1
dlaln2.f
Go to the documentation of this file.
1  SUBROUTINE dlaln2( LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B,
2  $ LDB, WR, WI, X, LDX, SCALE, XNORM, INFO )
3 *
4 * -- LAPACK auxiliary routine (version 3.0) --
5 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
6 * Courant Institute, Argonne National Lab, and Rice University
7 * October 31, 1992
8 *
9 * .. Scalar Arguments ..
10  LOGICAL LTRANS
11  INTEGER INFO, LDA, LDB, LDX, NA, NW
12  DOUBLE PRECISION CA, D1, D2, SCALE, SMIN, WI, WR, XNORM
13 * ..
14 * .. Array Arguments ..
15  DOUBLE PRECISION A( LDA, * ), B( LDB, * ), X( LDX, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DLALN2 solves a system of the form (ca A - w D ) X = s B
22 * or (ca A' - w D) X = s B with possible scaling ("s") and
23 * perturbation of A. (A' means A-transpose.)
24 *
25 * A is an NA x NA real matrix, ca is a real scalar, D is an NA x NA
26 * real diagonal matrix, w is a real or complex value, and X and B are
27 * NA x 1 matrices -- real if w is real, complex if w is complex. NA
28 * may be 1 or 2.
29 *
30 * If w is complex, X and B are represented as NA x 2 matrices,
31 * the first column of each being the real part and the second
32 * being the imaginary part.
33 *
34 * "s" is a scaling factor (.LE. 1), computed by DLALN2, which is
35 * so chosen that X can be computed without overflow. X is further
36 * scaled if necessary to assure that norm(ca A - w D)*norm(X) is less
37 * than overflow.
38 *
39 * If both singular values of (ca A - w D) are less than SMIN,
40 * SMIN*identity will be used instead of (ca A - w D). If only one
41 * singular value is less than SMIN, one element of (ca A - w D) will be
42 * perturbed enough to make the smallest singular value roughly SMIN.
43 * If both singular values are at least SMIN, (ca A - w D) will not be
44 * perturbed. In any case, the perturbation will be at most some small
45 * multiple of max( SMIN, ulp*norm(ca A - w D) ). The singular values
46 * are computed by infinity-norm approximations, and thus will only be
47 * correct to a factor of 2 or so.
48 *
49 * Note: all input quantities are assumed to be smaller than overflow
50 * by a reasonable factor. (See BIGNUM.)
51 *
52 * Arguments
53 * ==========
54 *
55 * LTRANS (input) LOGICAL
56 * =.TRUE.: A-transpose will be used.
57 * =.FALSE.: A will be used (not transposed.)
58 *
59 * NA (input) INTEGER
60 * The size of the matrix A. It may (only) be 1 or 2.
61 *
62 * NW (input) INTEGER
63 * 1 if "w" is real, 2 if "w" is complex. It may only be 1
64 * or 2.
65 *
66 * SMIN (input) DOUBLE PRECISION
67 * The desired lower bound on the singular values of A. This
68 * should be a safe distance away from underflow or overflow,
69 * say, between (underflow/machine precision) and (machine
70 * precision * overflow ). (See BIGNUM and ULP.)
71 *
72 * CA (input) DOUBLE PRECISION
73 * The coefficient c, which A is multiplied by.
74 *
75 * A (input) DOUBLE PRECISION array, dimension (LDA,NA)
76 * The NA x NA matrix A.
77 *
78 * LDA (input) INTEGER
79 * The leading dimension of A. It must be at least NA.
80 *
81 * D1 (input) DOUBLE PRECISION
82 * The 1,1 element in the diagonal matrix D.
83 *
84 * D2 (input) DOUBLE PRECISION
85 * The 2,2 element in the diagonal matrix D. Not used if NW=1.
86 *
87 * B (input) DOUBLE PRECISION array, dimension (LDB,NW)
88 * The NA x NW matrix B (right-hand side). If NW=2 ("w" is
89 * complex), column 1 contains the real part of B and column 2
90 * contains the imaginary part.
91 *
92 * LDB (input) INTEGER
93 * The leading dimension of B. It must be at least NA.
94 *
95 * WR (input) DOUBLE PRECISION
96 * The real part of the scalar "w".
97 *
98 * WI (input) DOUBLE PRECISION
99 * 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
dladiv
subroutine dladiv(A, B, C, D, P, Q)
Definition: dladiv.f:2
dlaln2
subroutine dlaln2(LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B, LDB, WR, WI, X, LDX, SCALE, XNORM, INFO)
Definition: dlaln2.f:3