1: #region Translated by Jose Antonio De Santiago-Castillo.
2:
3: //Translated by Jose Antonio De Santiago-Castillo.
4: //E-mail:JAntonioDeSantiago@gmail.com
5: //Web: www.DotNumerics.com
6: //
7: //Fortran to C# Translation.
8: //Translated by:
9: //F2CSharp Version 0.71 (November 10, 2009)
10: //Code Optimizations: None
11: //
12: #endregion
13:
14: using System;
15: using DotNumerics.FortranLibrary;
16:
17: namespace DotNumerics.CSLapack
18: {
19: /// <summary>
20: /// -- LAPACK auxiliary routine (version 3.1) --
21: /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
22: /// November 2006
23: /// Purpose
24: /// =======
25: ///
26: /// DLALN2 solves a system of the form (ca A - w D ) X = s B
27: /// or (ca A' - w D) X = s B with possible scaling ("s") and
28: /// perturbation of A. (A' means A-transpose.)
29: ///
30: /// A is an NA x NA real matrix, ca is a real scalar, D is an NA x NA
31: /// real diagonal matrix, w is a real or complex value, and X and B are
32: /// NA x 1 matrices -- real if w is real, complex if w is complex. NA
33: /// may be 1 or 2.
34: ///
35: /// If w is complex, X and B are represented as NA x 2 matrices,
36: /// the first column of each being the real part and the second
37: /// being the imaginary part.
38: ///
39: /// "s" is a scaling factor (.LE. 1), computed by DLALN2, which is
40: /// so chosen that X can be computed without overflow. X is further
41: /// scaled if necessary to assure that norm(ca A - w D)*norm(X) is less
42: /// than overflow.
43: ///
44: /// If both singular values of (ca A - w D) are less than SMIN,
45: /// SMIN*identity will be used instead of (ca A - w D). If only one
46: /// singular value is less than SMIN, one element of (ca A - w D) will be
47: /// perturbed enough to make the smallest singular value roughly SMIN.
48: /// If both singular values are at least SMIN, (ca A - w D) will not be
49: /// perturbed. In any case, the perturbation will be at most some small
50: /// multiple of max( SMIN, ulp*norm(ca A - w D) ). The singular values
51: /// are computed by infinity-norm approximations, and thus will only be
52: /// correct to a factor of 2 or so.
53: ///
54: /// Note: all input quantities are assumed to be smaller than overflow
55: /// by a reasonable factor. (See BIGNUM.)
56: ///
57: ///</summary>
58: public class DLALN2
59: {
60:
61:
62: #region Dependencies
63:
64: DLAMCH _dlamch; DLADIV _dladiv;
65:
66: #endregion
67:
68:
69: #region Fields
70:
71: const double ZERO = 0.0E0; const double ONE = 1.0E0; const double TWO = 2.0E0; int ICMAX = 0; int J = 0; double BBND = 0;
72: double BI1 = 0;double BI2 = 0; double BIGNUM = 0; double BNORM = 0; double BR1 = 0; double BR2 = 0; double CI21 = 0;
73: double CI22 = 0;double CMAX = 0; double CNORM = 0; double CR21 = 0; double CR22 = 0; double CSI = 0; double CSR = 0;
74: double LI21 = 0;double LR21 = 0; double SMINI = 0; double SMLNUM = 0; double TEMP = 0; double U22ABS = 0; double UI11 = 0;
75: double UI11R = 0;double UI12 = 0; double UI12S = 0; double UI22 = 0; double UR11 = 0; double UR11R = 0; double UR12 = 0;
76: double UR12S = 0;double UR22 = 0; double XI1 = 0; double XI2 = 0; double XR1 = 0; double XR2 = 0;
77: bool[] RSWAP = new bool[4]; int o_rswap = -1;bool[] ZSWAP = new bool[4]; int o_zswap = -1;
78: int[] IPIVOT = new int[4 * 4]; int o_ipivot = -5;double[] CIV = new double[4]; int o_civ = -1;
79: double[] CRV = new double[4]; int o_crv = -1;
80:
81: #endregion
82:
83: public DLALN2(DLAMCH dlamch, DLADIV dladiv)
84: {
85:
86:
87: #region Set Dependencies
88:
89: this._dlamch = dlamch; this._dladiv = dladiv;
90:
91: #endregion
92:
93:
94: #region Data Inicializacion
95:
96: //ZSWAP/.FALSE.,.FALSE.,.TRUE.,.TRUE.
97: ZSWAP[1 + o_zswap] = false;
98: ZSWAP[2 + o_zswap] = false;
99: ZSWAP[3 + o_zswap] = true;
100: ZSWAP[4 + o_zswap] = true;
101: //RSWAP/.FALSE.,.TRUE.,.FALSE.,.TRUE.
102: RSWAP[1 + o_rswap] = false;
103: RSWAP[2 + o_rswap] = true;
104: RSWAP[3 + o_rswap] = false;
105: RSWAP[4 + o_rswap] = true;
106: //IPIVOT/1,2,3,4,2,1,4,3,3,4,1,2,4,3,2,1
107: IPIVOT[0] = 1;
108: IPIVOT[1] = 2;
109: IPIVOT[2] = 3;
110: IPIVOT[3] = 4;
111: IPIVOT[4] = 2;
112: IPIVOT[5] = 1;
113: IPIVOT[6] = 4;
114: IPIVOT[7] = 3;
115: IPIVOT[8] = 3;
116: IPIVOT[9] = 4;
117: IPIVOT[10] = 1;
118: IPIVOT[11] = 2;
119: IPIVOT[12] = 4;
120: IPIVOT[12] = 3;
121: IPIVOT[14] = 2;
122: IPIVOT[15] = 1;
123:
124: #endregion
125:
126: }
127:
128: public DLALN2()
129: {
130:
131:
132: #region Dependencies (Initialization)
133:
134: LSAME lsame = new LSAME();
135: DLAMC3 dlamc3 = new DLAMC3();
136: DLADIV dladiv = new DLADIV();
137: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
138: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
139: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
140: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
141: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
142:
143: #endregion
144:
145:
146: #region Set Dependencies
147:
148: this._dlamch = dlamch; this._dladiv = dladiv;
149:
150: #endregion
151:
152:
153: #region Data Inicializacion
154:
155: //ZSWAP/.FALSE.,.FALSE.,.TRUE.,.TRUE.
156: ZSWAP[1 + o_zswap] = false;
157: ZSWAP[2 + o_zswap] = false;
158: ZSWAP[3 + o_zswap] = true;
159: ZSWAP[4 + o_zswap] = true;
160: //RSWAP/.FALSE.,.TRUE.,.FALSE.,.TRUE.
161: RSWAP[1 + o_rswap] = false;
162: RSWAP[2 + o_rswap] = true;
163: RSWAP[3 + o_rswap] = false;
164: RSWAP[4 + o_rswap] = true;
165: //IPIVOT/1,2,3,4,2,1,4,3,3,4,1,2,4,3,2,1
166: IPIVOT[0] = 1;
167: IPIVOT[1] = 2;
168: IPIVOT[2] = 3;
169: IPIVOT[3] = 4;
170: IPIVOT[4] = 2;
171: IPIVOT[5] = 1;
172: IPIVOT[6] = 4;
173: IPIVOT[7] = 3;
174: IPIVOT[8] = 3;
175: IPIVOT[9] = 4;
176: IPIVOT[10] = 1;
177: IPIVOT[11] = 2;
178: IPIVOT[12] = 4;
179: IPIVOT[12] = 3;
180: IPIVOT[14] = 2;
181: IPIVOT[15] = 1;
182:
183: #endregion
184:
185: }
186: /// <summary>
187: /// Purpose
188: /// =======
189: ///
190: /// DLALN2 solves a system of the form (ca A - w D ) X = s B
191: /// or (ca A' - w D) X = s B with possible scaling ("s") and
192: /// perturbation of A. (A' means A-transpose.)
193: ///
194: /// A is an NA x NA real matrix, ca is a real scalar, D is an NA x NA
195: /// real diagonal matrix, w is a real or complex value, and X and B are
196: /// NA x 1 matrices -- real if w is real, complex if w is complex. NA
197: /// may be 1 or 2.
198: ///
199: /// If w is complex, X and B are represented as NA x 2 matrices,
200: /// the first column of each being the real part and the second
201: /// being the imaginary part.
202: ///
203: /// "s" is a scaling factor (.LE. 1), computed by DLALN2, which is
204: /// so chosen that X can be computed without overflow. X is further
205: /// scaled if necessary to assure that norm(ca A - w D)*norm(X) is less
206: /// than overflow.
207: ///
208: /// If both singular values of (ca A - w D) are less than SMIN,
209: /// SMIN*identity will be used instead of (ca A - w D). If only one
210: /// singular value is less than SMIN, one element of (ca A - w D) will be
211: /// perturbed enough to make the smallest singular value roughly SMIN.
212: /// If both singular values are at least SMIN, (ca A - w D) will not be
213: /// perturbed. In any case, the perturbation will be at most some small
214: /// multiple of max( SMIN, ulp*norm(ca A - w D) ). The singular values
215: /// are computed by infinity-norm approximations, and thus will only be
216: /// correct to a factor of 2 or so.
217: ///
218: /// Note: all input quantities are assumed to be smaller than overflow
219: /// by a reasonable factor. (See BIGNUM.)
220: ///
221: ///</summary>
222: /// <param name="LTRANS">
223: /// (input) LOGICAL
224: /// =.TRUE.: A-transpose will be used.
225: /// =.FALSE.: A will be used (not transposed.)
226: ///</param>
227: /// <param name="NA">
228: /// (input) INTEGER
229: /// The size of the matrix A. It may (only) be 1 or 2.
230: ///</param>
231: /// <param name="NW">
232: /// (input) INTEGER
233: /// 1 if "w" is real, 2 if "w" is complex. It may only be 1
234: /// or 2.
235: ///</param>
236: /// <param name="SMIN">
237: /// (input) DOUBLE PRECISION
238: /// The desired lower bound on the singular values of A. This
239: /// should be a safe distance away from underflow or overflow,
240: /// say, between (underflow/machine precision) and (machine
241: /// precision * overflow ). (See BIGNUM and ULP.)
242: ///</param>
243: /// <param name="CA">
244: /// (input) DOUBLE PRECISION
245: /// The coefficient c, which A is multiplied by.
246: ///</param>
247: /// <param name="A">
248: /// is an NA x NA real matrix, ca is a real scalar, D is an NA x NA
249: ///</param>
250: /// <param name="LDA">
251: /// (input) INTEGER
252: /// The leading dimension of A. It must be at least NA.
253: ///</param>
254: /// <param name="D1">
255: /// (input) DOUBLE PRECISION
256: /// The 1,1 element in the diagonal matrix D.
257: ///</param>
258: /// <param name="D2">
259: /// (input) DOUBLE PRECISION
260: /// The 2,2 element in the diagonal matrix D. Not used if NW=1.
261: ///</param>
262: /// <param name="B">
263: /// (input) DOUBLE PRECISION array, dimension (LDB,NW)
264: /// The NA x NW matrix B (right-hand side). If NW=2 ("w" is
265: /// complex), column 1 contains the real part of B and column 2
266: /// contains the imaginary part.
267: ///</param>
268: /// <param name="LDB">
269: /// (input) INTEGER
270: /// The leading dimension of B. It must be at least NA.
271: ///</param>
272: /// <param name="WR">
273: /// (input) DOUBLE PRECISION
274: /// The real part of the scalar "w".
275: ///</param>
276: /// <param name="WI">
277: /// (input) DOUBLE PRECISION
278: /// The imaginary part of the scalar "w". Not used if NW=1.
279: ///</param>
280: /// <param name="X">
281: /// (output) DOUBLE PRECISION array, dimension (LDX,NW)
282: /// The NA x NW matrix X (unknowns), as computed by DLALN2.
283: /// If NW=2 ("w" is complex), on exit, column 1 will contain
284: /// the real part of X and column 2 will contain the imaginary
285: /// part.
286: ///</param>
287: /// <param name="LDX">
288: /// (input) INTEGER
289: /// The leading dimension of X. It must be at least NA.
290: ///</param>
291: /// <param name="SCALE">
292: /// (output) DOUBLE PRECISION
293: /// The scale factor that B must be multiplied by to insure
294: /// that overflow does not occur when computing X. Thus,
295: /// (ca A - w D) X will be SCALE*B, not B (ignoring
296: /// perturbations of A.) It will be at most 1.
297: ///</param>
298: /// <param name="XNORM">
299: /// (output) DOUBLE PRECISION
300: /// The infinity-norm of X, when X is regarded as an NA x NW
301: /// real matrix.
302: ///</param>
303: /// <param name="INFO">
304: /// (output) INTEGER
305: /// An error flag. It will be set to zero if no error occurs,
306: /// a negative number if an argument is in error, or a positive
307: /// number if ca A - w D had to be perturbed.
308: /// The possible values are:
309: /// = 0: No error occurred, and (ca A - w D) did not have to be
310: /// perturbed.
311: /// = 1: (ca A - w D) had to be perturbed to make its smallest
312: /// (or only) singular value greater than SMIN.
313: /// NOTE: In the interests of speed, this routine does not
314: /// check the inputs for errors.
315: ///</param>
316: public void Run(bool LTRANS, int NA, int NW, double SMIN, double CA, double[] A, int offset_a
317: , int LDA, double D1, double D2, double[] B, int offset_b, int LDB, double WR
318: , double WI, ref double[] X, int offset_x, int LDX, ref double SCALE, ref double XNORM, ref int INFO)
319: {
320:
321: #region Array Index Correction
322:
323: int o_a = -1 - LDA + offset_a; int o_b = -1 - LDB + offset_b; int o_x = -1 - LDX + offset_x;
324:
325: #endregion
326:
327:
328: #region Prolog
329:
330:
331:
332: // *
333: // * SUBRUTINA MODIFICADA, SE ELIMINO LA EQUIVALENCIA
334: // *
335: // * EQUIVALENCE ( CI( 1, 1 ), CIV( 1 ) ),
336: // * $ ( CR( 1, 1 ), CRV( 1 ) )
337: // *
338: // * DEBIDO A QUE EL TRADUCTOR NO LA SOPORTA
339:
340: // *
341: // * -- LAPACK auxiliary routine (version 3.1) --
342: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
343: // * November 2006
344: // *
345: // * .. Scalar Arguments ..
346: // * ..
347: // * .. Array Arguments ..
348: // * ..
349: // *
350: // * Purpose
351: // * =======
352: // *
353: // * DLALN2 solves a system of the form (ca A - w D ) X = s B
354: // * or (ca A' - w D) X = s B with possible scaling ("s") and
355: // * perturbation of A. (A' means A-transpose.)
356: // *
357: // * A is an NA x NA real matrix, ca is a real scalar, D is an NA x NA
358: // * real diagonal matrix, w is a real or complex value, and X and B are
359: // * NA x 1 matrices -- real if w is real, complex if w is complex. NA
360: // * may be 1 or 2.
361: // *
362: // * If w is complex, X and B are represented as NA x 2 matrices,
363: // * the first column of each being the real part and the second
364: // * being the imaginary part.
365: // *
366: // * "s" is a scaling factor (.LE. 1), computed by DLALN2, which is
367: // * so chosen that X can be computed without overflow. X is further
368: // * scaled if necessary to assure that norm(ca A - w D)*norm(X) is less
369: // * than overflow.
370: // *
371: // * If both singular values of (ca A - w D) are less than SMIN,
372: // * SMIN*identity will be used instead of (ca A - w D). If only one
373: // * singular value is less than SMIN, one element of (ca A - w D) will be
374: // * perturbed enough to make the smallest singular value roughly SMIN.
375: // * If both singular values are at least SMIN, (ca A - w D) will not be
376: // * perturbed. In any case, the perturbation will be at most some small
377: // * multiple of max( SMIN, ulp*norm(ca A - w D) ). The singular values
378: // * are computed by infinity-norm approximations, and thus will only be
379: // * correct to a factor of 2 or so.
380: // *
381: // * Note: all input quantities are assumed to be smaller than overflow
382: // * by a reasonable factor. (See BIGNUM.)
383: // *
384: // * Arguments
385: // * ==========
386: // *
387: // * LTRANS (input) LOGICAL
388: // * =.TRUE.: A-transpose will be used.
389: // * =.FALSE.: A will be used (not transposed.)
390: // *
391: // * NA (input) INTEGER
392: // * The size of the matrix A. It may (only) be 1 or 2.
393: // *
394: // * NW (input) INTEGER
395: // * 1 if "w" is real, 2 if "w" is complex. It may only be 1
396: // * or 2.
397: // *
398: // * SMIN (input) DOUBLE PRECISION
399: // * The desired lower bound on the singular values of A. This
400: // * should be a safe distance away from underflow or overflow,
401: // * say, between (underflow/machine precision) and (machine
402: // * precision * overflow ). (See BIGNUM and ULP.)
403: // *
404: // * CA (input) DOUBLE PRECISION
405: // * The coefficient c, which A is multiplied by.
406: // *
407: // * A (input) DOUBLE PRECISION array, dimension (LDA,NA)
408: // * The NA x NA matrix A.
409: // *
410: // * LDA (input) INTEGER
411: // * The leading dimension of A. It must be at least NA.
412: // *
413: // * D1 (input) DOUBLE PRECISION
414: // * The 1,1 element in the diagonal matrix D.
415: // *
416: // * D2 (input) DOUBLE PRECISION
417: // * The 2,2 element in the diagonal matrix D. Not used if NW=1.
418: // *
419: // * B (input) DOUBLE PRECISION array, dimension (LDB,NW)
420: // * The NA x NW matrix B (right-hand side). If NW=2 ("w" is
421: // * complex), column 1 contains the real part of B and column 2
422: // * contains the imaginary part.
423: // *
424: // * LDB (input) INTEGER
425: // * The leading dimension of B. It must be at least NA.
426: // *
427: // * WR (input) DOUBLE PRECISION
428: // * The real part of the scalar "w".
429: // *
430: // * WI (input) DOUBLE PRECISION
431: // * The imaginary part of the scalar "w". Not used if NW=1.
432: // *
433: // * X (output) DOUBLE PRECISION array, dimension (LDX,NW)
434: // * The NA x NW matrix X (unknowns), as computed by DLALN2.
435: // * If NW=2 ("w" is complex), on exit, column 1 will contain
436: // * the real part of X and column 2 will contain the imaginary
437: // * part.
438: // *
439: // * LDX (input) INTEGER
440: // * The leading dimension of X. It must be at least NA.
441: // *
442: // * SCALE (output) DOUBLE PRECISION
443: // * The scale factor that B must be multiplied by to insure
444: // * that overflow does not occur when computing X. Thus,
445: // * (ca A - w D) X will be SCALE*B, not B (ignoring
446: // * perturbations of A.) It will be at most 1.
447: // *
448: // * XNORM (output) DOUBLE PRECISION
449: // * The infinity-norm of X, when X is regarded as an NA x NW
450: // * real matrix.
451: // *
452: // * INFO (output) INTEGER
453: // * An error flag. It will be set to zero if no error occurs,
454: // * a negative number if an argument is in error, or a positive
455: // * number if ca A - w D had to be perturbed.
456: // * The possible values are:
457: // * = 0: No error occurred, and (ca A - w D) did not have to be
458: // * perturbed.
459: // * = 1: (ca A - w D) had to be perturbed to make its smallest
460: // * (or only) singular value greater than SMIN.
461: // * NOTE: In the interests of speed, this routine does not
462: // * check the inputs for errors.
463: // *
464: // * =====================================================================
465: // *
466: // * .. Parameters ..
467: // * ..
468: // * .. Local Scalars ..
469: // * ..
470: // * .. Local Arrays ..
471:
472:
473:
474: // *-----------------------------INICIA MI MODIFICACION--------------------------
475:
476: // * DOUBLE PRECISION CI( 2, 2 ), CIV( 4 ), CR( 2, 2 ), CRV( 4 )
477: // *-----------------------------TERMINA MI MODIFICACION-------------------------
478:
479:
480:
481: // * ..
482: // * .. External Functions ..
483: // * ..
484: // * .. External Subroutines ..
485: // * ..
486: // * .. Intrinsic Functions ..
487: // INTRINSIC ABS, MAX;
488: // * ..
489: // * .. Equivalences ..
490:
491:
492: // *-----------------------------INICIA MI MODIFICACION--------------------------
493: // * EQUIVALENCE ( CI( 1, 1 ), CIV( 1 ) ),
494: // * $ ( CR( 1, 1 ), CRV( 1 ) )
495: // *-----------------------------TERMINA MI MODIFICACION-------------------------
496:
497:
498:
499: // * .. Data statements ..
500: // * ..
501: // * .. Executable Statements ..
502: // *
503: // * Compute BIGNUM
504: // *
505:
506: #endregion
507:
508:
509: #region Body
510:
511: SMLNUM = TWO * this._dlamch.Run("Safe minimum");
512: BIGNUM = ONE / SMLNUM;
513: SMINI = Math.Max(SMIN, SMLNUM);
514: // *
515: // * Don't check for input errors
516: // *
517: INFO = 0;
518: // *
519: // * Standard Initializations
520: // *
521: SCALE = ONE;
522: // *
523: if (NA == 1)
524: {
525: // *
526: // * 1 x 1 (i.e., scalar) system C X = B
527: // *
528: if (NW == 1)
529: {
530: // *
531: // * Real 1x1 system.
532: // *
533: // * C = ca A - w D
534: // *
535: CSR = CA * A[1+1 * LDA + o_a] - WR * D1;
536: CNORM = Math.Abs(CSR);
537: // *
538: // * If | C | < SMINI, use C = SMINI
539: // *
540: if (CNORM < SMINI)
541: {
542: CSR = SMINI;
543: CNORM = SMINI;
544: INFO = 1;
545: }
546: // *
547: // * Check scaling for X = B / C
548: // *
549: BNORM = Math.Abs(B[1+1 * LDB + o_b]);
550: if (CNORM < ONE && BNORM > ONE)
551: {
552: if (BNORM > BIGNUM * CNORM) SCALE = ONE / BNORM;
553: }
554: // *
555: // * Compute X
556: // *
557: X[1+1 * LDX + o_x] = (B[1+1 * LDB + o_b] * SCALE) / CSR;
558: XNORM = Math.Abs(X[1+1 * LDX + o_x]);
559: }
560: else
561: {
562: // *
563: // * Complex 1x1 system (w is complex)
564: // *
565: // * C = ca A - w D
566: // *
567: CSR = CA * A[1+1 * LDA + o_a] - WR * D1;
568: CSI = - WI * D1;
569: CNORM = Math.Abs(CSR) + Math.Abs(CSI);
570: // *
571: // * If | C | < SMINI, use C = SMINI
572: // *
573: if (CNORM < SMINI)
574: {
575: CSR = SMINI;
576: CSI = ZERO;
577: CNORM = SMINI;
578: INFO = 1;
579: }
580: // *
581: // * Check scaling for X = B / C
582: // *
583: BNORM = Math.Abs(B[1+1 * LDB + o_b]) + Math.Abs(B[1+2 * LDB + o_b]);
584: if (CNORM < ONE && BNORM > ONE)
585: {
586: if (BNORM > BIGNUM * CNORM) SCALE = ONE / BNORM;
587: }
588: // *
589: // * Compute X
590: // *
591: this._dladiv.Run(SCALE * B[1+1 * LDB + o_b], SCALE * B[1+2 * LDB + o_b], CSR, CSI, ref X[1+1 * LDX + o_x], ref X[1+2 * LDX + o_x]);
592: XNORM = Math.Abs(X[1+1 * LDX + o_x]) + Math.Abs(X[1+2 * LDX + o_x]);
593: }
594: // *
595: }
596: else
597: {
598: // *
599: // * 2x2 System
600: // *
601: // * Compute the real part of C = ca A - w D (or ca A' - w D )
602: // *
603:
604:
605:
606:
607: // *-----------------------------INICIA MI MODIFICACION--------------------------
608:
609: // * CR( 1, 1 ) = CA*A( 1, 1 ) - WR*D1
610: // * CR( 2, 2 ) = CA*A( 2, 2 ) - WR*D2
611:
612: CRV[1 + o_crv] = CA * A[1+1 * LDA + o_a] - WR * D1;
613: CRV[4 + o_crv] = CA * A[2+2 * LDA + o_a] - WR * D2;
614:
615: if (LTRANS)
616: {
617: // * CR( 1, 2 ) = CA*A( 2, 1 )
618: // * CR( 2, 1 ) = CA*A( 1, 2 )
619:
620: CRV[3 + o_crv] = CA * A[2+1 * LDA + o_a];
621: CRV[2 + o_crv] = CA * A[1+2 * LDA + o_a];
622:
623: }
624: else
625: {
626: // * CR( 2, 1 ) = CA*A( 2, 1 )
627: // * CR( 1, 2 ) = CA*A( 1, 2 )
628:
629:
630: CRV[2 + o_crv] = CA * A[2+1 * LDA + o_a];
631: CRV[3 + o_crv] = CA * A[1+2 * LDA + o_a];
632:
633: }
634:
635: // *-----------------------------TERMINA MI MODIFICACION-------------------------
636:
637:
638:
639: // *
640: if (NW == 1)
641: {
642: // *
643: // * Real 2x2 system (w is real)
644: // *
645: // * Find the largest element in C
646: // *
647: CMAX = ZERO;
648: ICMAX = 0;
649: // *
650: for (J = 1; J <= 4; J++)
651: {
652: if (Math.Abs(CRV[J + o_crv]) > CMAX)
653: {
654: CMAX = Math.Abs(CRV[J + o_crv]);
655: ICMAX = J;
656: }
657: }
658: // *
659: // * If norm(C) < SMINI, use SMINI*identity.
660: // *
661: if (CMAX < SMINI)
662: {
663: BNORM = Math.Max(Math.Abs(B[1+1 * LDB + o_b]), Math.Abs(B[2+1 * LDB + o_b]));
664: if (SMINI < ONE && BNORM > ONE)
665: {
666: if (BNORM > BIGNUM * SMINI) SCALE = ONE / BNORM;
667: }
668: TEMP = SCALE / SMINI;
669: X[1+1 * LDX + o_x] = TEMP * B[1+1 * LDB + o_b];
670: X[2+1 * LDX + o_x] = TEMP * B[2+1 * LDB + o_b];
671: XNORM = TEMP * BNORM;
672: INFO = 1;
673: return;
674: }
675: // *
676: // * Gaussian elimination with complete pivoting.
677: // *
678: UR11 = CRV[ICMAX + o_crv];
679: CR21 = CRV[IPIVOT[2+ICMAX * 4 + o_ipivot] + o_crv];
680: UR12 = CRV[IPIVOT[3+ICMAX * 4 + o_ipivot] + o_crv];
681: CR22 = CRV[IPIVOT[4+ICMAX * 4 + o_ipivot] + o_crv];
682: UR11R = ONE / UR11;
683: LR21 = UR11R * CR21;
684: UR22 = CR22 - UR12 * LR21;
685: // *
686: // * If smaller pivot < SMINI, use SMINI
687: // *
688: if (Math.Abs(UR22) < SMINI)
689: {
690: UR22 = SMINI;
691: INFO = 1;
692: }
693: if (RSWAP[ICMAX + o_rswap])
694: {
695: BR1 = B[2+1 * LDB + o_b];
696: BR2 = B[1+1 * LDB + o_b];
697: }
698: else
699: {
700: BR1 = B[1+1 * LDB + o_b];
701: BR2 = B[2+1 * LDB + o_b];
702: }
703: BR2 = BR2 - LR21 * BR1;
704: BBND = Math.Max(Math.Abs(BR1 * (UR22 * UR11R)), Math.Abs(BR2));
705: if (BBND > ONE && Math.Abs(UR22) < ONE)
706: {
707: if (BBND >= BIGNUM * Math.Abs(UR22)) SCALE = ONE / BBND;
708: }
709: // *
710: XR2 = (BR2 * SCALE) / UR22;
711: XR1 = (SCALE * BR1) * UR11R - XR2 * (UR11R * UR12);
712: if (ZSWAP[ICMAX + o_zswap])
713: {
714: X[1+1 * LDX + o_x] = XR2;
715: X[2+1 * LDX + o_x] = XR1;
716: }
717: else
718: {
719: X[1+1 * LDX + o_x] = XR1;
720: X[2+1 * LDX + o_x] = XR2;
721: }
722: XNORM = Math.Max(Math.Abs(XR1), Math.Abs(XR2));
723: // *
724: // * Further scaling if norm(A) norm(X) > overflow
725: // *
726: if (XNORM > ONE && CMAX > ONE)
727: {
728: if (XNORM > BIGNUM / CMAX)
729: {
730: TEMP = CMAX / BIGNUM;
731: X[1+1 * LDX + o_x] = TEMP * X[1+1 * LDX + o_x];
732: X[2+1 * LDX + o_x] = TEMP * X[2+1 * LDX + o_x];
733: XNORM = TEMP * XNORM;
734: SCALE = TEMP * SCALE;
735: }
736: }
737: }
738: else
739: {
740: // *
741: // * Complex 2x2 system (w is complex)
742: // *
743: // * Find the largest element in C
744: // *
745:
746:
747: // *-----------------------------INICIA MI MODIFICACION--------------------------
748: // * CI( 1, 1 ) = -WI*D1
749: // * CI( 2, 1 ) = ZERO
750: // * CI( 1, 2 ) = ZERO
751: // * CI( 2, 2 ) = -WI*D2
752: CIV[1 + o_civ] = - WI * D1;
753: CIV[2 + o_civ] = ZERO;
754: CIV[3 + o_civ] = ZERO;
755: CIV[4 + o_civ] = - WI * D2;
756: // *-----------------------------TERMINA MI MODIFICACION-------------------------
757:
758:
759:
760:
761:
762: CMAX = ZERO;
763: ICMAX = 0;
764: // *
765: for (J = 1; J <= 4; J++)
766: {
767: if (Math.Abs(CRV[J + o_crv]) + Math.Abs(CIV[J + o_civ]) > CMAX)
768: {
769: CMAX = Math.Abs(CRV[J + o_crv]) + Math.Abs(CIV[J + o_civ]);
770: ICMAX = J;
771: }
772: }
773: // *
774: // * If norm(C) < SMINI, use SMINI*identity.
775: // *
776: if (CMAX < SMINI)
777: {
778: BNORM = Math.Max(Math.Abs(B[1+1 * LDB + o_b]) + Math.Abs(B[1+2 * LDB + o_b]), Math.Abs(B[2+1 * LDB + o_b]) + Math.Abs(B[2+2 * LDB + o_b]));
779: if (SMINI < ONE && BNORM > ONE)
780: {
781: if (BNORM > BIGNUM * SMINI) SCALE = ONE / BNORM;
782: }
783: TEMP = SCALE / SMINI;
784: X[1+1 * LDX + o_x] = TEMP * B[1+1 * LDB + o_b];
785: X[2+1 * LDX + o_x] = TEMP * B[2+1 * LDB + o_b];
786: X[1+2 * LDX + o_x] = TEMP * B[1+2 * LDB + o_b];
787: X[2+2 * LDX + o_x] = TEMP * B[2+2 * LDB + o_b];
788: XNORM = TEMP * BNORM;
789: INFO = 1;
790: return;
791: }
792: // *
793: // * Gaussian elimination with complete pivoting.
794: // *
795: UR11 = CRV[ICMAX + o_crv];
796: UI11 = CIV[ICMAX + o_civ];
797: CR21 = CRV[IPIVOT[2+ICMAX * 4 + o_ipivot] + o_crv];
798: CI21 = CIV[IPIVOT[2+ICMAX * 4 + o_ipivot] + o_civ];
799: UR12 = CRV[IPIVOT[3+ICMAX * 4 + o_ipivot] + o_crv];
800: UI12 = CIV[IPIVOT[3+ICMAX * 4 + o_ipivot] + o_civ];
801: CR22 = CRV[IPIVOT[4+ICMAX * 4 + o_ipivot] + o_crv];
802: CI22 = CIV[IPIVOT[4+ICMAX * 4 + o_ipivot] + o_civ];
803: if (ICMAX == 1 || ICMAX == 4)
804: {
805: // *
806: // * Code when off-diagonals of pivoted C are real
807: // *
808: if (Math.Abs(UR11) > Math.Abs(UI11))
809: {
810: TEMP = UI11 / UR11;
811: UR11R = ONE / (UR11 * (ONE + Math.Pow(TEMP,2)));
812: UI11R = - TEMP * UR11R;
813: }
814: else
815: {
816: TEMP = UR11 / UI11;
817: UI11R = - ONE / (UI11 * (ONE + Math.Pow(TEMP,2)));
818: UR11R = - TEMP * UI11R;
819: }
820: LR21 = CR21 * UR11R;
821: LI21 = CR21 * UI11R;
822: UR12S = UR12 * UR11R;
823: UI12S = UR12 * UI11R;
824: UR22 = CR22 - UR12 * LR21;
825: UI22 = CI22 - UR12 * LI21;
826: }
827: else
828: {
829: // *
830: // * Code when diagonals of pivoted C are real
831: // *
832: UR11R = ONE / UR11;
833: UI11R = ZERO;
834: LR21 = CR21 * UR11R;
835: LI21 = CI21 * UR11R;
836: UR12S = UR12 * UR11R;
837: UI12S = UI12 * UR11R;
838: UR22 = CR22 - UR12 * LR21 + UI12 * LI21;
839: UI22 = - UR12 * LI21 - UI12 * LR21;
840: }
841: U22ABS = Math.Abs(UR22) + Math.Abs(UI22);
842: // *
843: // * If smaller pivot < SMINI, use SMINI
844: // *
845: if (U22ABS < SMINI)
846: {
847: UR22 = SMINI;
848: UI22 = ZERO;
849: INFO = 1;
850: }
851: if (RSWAP[ICMAX + o_rswap])
852: {
853: BR2 = B[1+1 * LDB + o_b];
854: BR1 = B[2+1 * LDB + o_b];
855: BI2 = B[1+2 * LDB + o_b];
856: BI1 = B[2+2 * LDB + o_b];
857: }
858: else
859: {
860: BR1 = B[1+1 * LDB + o_b];
861: BR2 = B[2+1 * LDB + o_b];
862: BI1 = B[1+2 * LDB + o_b];
863: BI2 = B[2+2 * LDB + o_b];
864: }
865: BR2 = BR2 - LR21 * BR1 + LI21 * BI1;
866: BI2 = BI2 - LI21 * BR1 - LR21 * BI1;
867: BBND = Math.Max((Math.Abs(BR1) + Math.Abs(BI1)) * (U22ABS * (Math.Abs(UR11R) + Math.Abs(UI11R))), Math.Abs(BR2) + Math.Abs(BI2));
868: if (BBND > ONE && U22ABS < ONE)
869: {
870: if (BBND >= BIGNUM * U22ABS)
871: {
872: SCALE = ONE / BBND;
873: BR1 = SCALE * BR1;
874: BI1 = SCALE * BI1;
875: BR2 = SCALE * BR2;
876: BI2 = SCALE * BI2;
877: }
878: }
879: // *
880: this._dladiv.Run(BR2, BI2, UR22, UI22, ref XR2, ref XI2);
881: XR1 = UR11R * BR1 - UI11R * BI1 - UR12S * XR2 + UI12S * XI2;
882: XI1 = UI11R * BR1 + UR11R * BI1 - UI12S * XR2 - UR12S * XI2;
883: if (ZSWAP[ICMAX + o_zswap])
884: {
885: X[1+1 * LDX + o_x] = XR2;
886: X[2+1 * LDX + o_x] = XR1;
887: X[1+2 * LDX + o_x] = XI2;
888: X[2+2 * LDX + o_x] = XI1;
889: }
890: else
891: {
892: X[1+1 * LDX + o_x] = XR1;
893: X[2+1 * LDX + o_x] = XR2;
894: X[1+2 * LDX + o_x] = XI1;
895: X[2+2 * LDX + o_x] = XI2;
896: }
897: XNORM = Math.Max(Math.Abs(XR1) + Math.Abs(XI1), Math.Abs(XR2) + Math.Abs(XI2));
898: // *
899: // * Further scaling if norm(A) norm(X) > overflow
900: // *
901: if (XNORM > ONE && CMAX > ONE)
902: {
903: if (XNORM > BIGNUM / CMAX)
904: {
905: TEMP = CMAX / BIGNUM;
906: X[1+1 * LDX + o_x] = TEMP * X[1+1 * LDX + o_x];
907: X[2+1 * LDX + o_x] = TEMP * X[2+1 * LDX + o_x];
908: X[1+2 * LDX + o_x] = TEMP * X[1+2 * LDX + o_x];
909: X[2+2 * LDX + o_x] = TEMP * X[2+2 * LDX + o_x];
910: XNORM = TEMP * XNORM;
911: SCALE = TEMP * SCALE;
912: }
913: }
914: }
915: }
916: // *
917: return;
918: // *
919: // * End of DLALN2
920: // *
921:
922: #endregion
923:
924: }
925: }
926: }