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: /// DLASY2 solves for the N1 by N2 matrix X, 1 .LE. N1,N2 .LE. 2, in
27: ///
28: /// op(TL)*X + ISGN*X*op(TR) = SCALE*B,
29: ///
30: /// where TL is N1 by N1, TR is N2 by N2, B is N1 by N2, and ISGN = 1 or
31: /// -1. op(T) = T or T', where T' denotes the transpose of T.
32: ///
33: ///</summary>
34: public class DLASY2
35: {
36:
37:
38: #region Dependencies
39:
40: IDAMAX _idamax; DLAMCH _dlamch; DCOPY _dcopy; DSWAP _dswap;
41:
42: #endregion
43:
44:
45: #region Fields
46:
47: const double ZERO = 0.0E+0; const double ONE = 1.0E+0; const double TWO = 2.0E+0; const double HALF = 0.5E+0;
48: const double EIGHT = 8.0E+0;bool BSWAP = false; bool XSWAP = false; int I = 0; int IP = 0; int IPIV = 0; int IPSV = 0;
49: int J = 0;int JP = 0; int JPSV = 0; int K = 0; double BET = 0; double EPS = 0; double GAM = 0; double L21 = 0;
50: double SGN = 0;double SMIN = 0; double SMLNUM = 0; double TAU1 = 0; double TEMP = 0; double U11 = 0; double U12 = 0;
51: double U22 = 0;double XMAX = 0; bool[] BSWPIV = new bool[4]; int o_bswpiv = -1;
52: bool[] XSWPIV = new bool[4]; int o_xswpiv = -1;int[] JPIV = new int[4]; int o_jpiv = -1;
53: int[] LOCL21 = new int[4]; int o_locl21 = -1;int[] LOCU12 = new int[4]; int o_locu12 = -1;
54: int[] LOCU22 = new int[4]; int o_locu22 = -1;double[] BTMP = new double[4]; int offset_btmp = 0; int o_btmp = -1;
55: double[] T16 = new double[4 * 4]; int offset_t16 = 0; int o_t16 = -5;double[] TMP = new double[4]; int offset_tmp = 0; int o_tmp = -1;
56: double[] X2 = new double[2]; int o_x2 = -1;
57:
58: #endregion
59:
60: public DLASY2(IDAMAX idamax, DLAMCH dlamch, DCOPY dcopy, DSWAP dswap)
61: {
62:
63:
64: #region Set Dependencies
65:
66: this._idamax = idamax; this._dlamch = dlamch; this._dcopy = dcopy; this._dswap = dswap;
67:
68: #endregion
69:
70:
71: #region Data Inicializacion
72:
73: //LOCU12/3,4,1,2
74: LOCU12[1 + o_locu12] = 3;
75: LOCU12[2 + o_locu12] = 4;
76: LOCU12[3 + o_locu12] = 1;
77: LOCU12[4 + o_locu12] = 2;
78: //LOCL21/2,1,4,3
79: LOCL21[1 + o_locl21] = 2;
80: LOCL21[2 + o_locl21] = 1;
81: LOCL21[3 + o_locl21] = 4;
82: LOCL21[4 + o_locl21] = 3;
83: //LOCU22/4,3,2,1
84: LOCU22[1 + o_locu22] = 4;
85: LOCU22[2 + o_locu22] = 3;
86: LOCU22[3 + o_locu22] = 2;
87: LOCU22[4 + o_locu22] = 1;
88: //XSWPIV/.FALSE.,.FALSE.,.TRUE.,.TRUE.
89: XSWPIV[1 + o_xswpiv] = false;
90: XSWPIV[2 + o_xswpiv] = false;
91: XSWPIV[3 + o_xswpiv] = true;
92: XSWPIV[4 + o_xswpiv] = true;
93: //BSWPIV/.FALSE.,.TRUE.,.FALSE.,.TRUE.
94: BSWPIV[1 + o_bswpiv] = false;
95: BSWPIV[2 + o_bswpiv] = true;
96: BSWPIV[3 + o_bswpiv] = false;
97: BSWPIV[4 + o_bswpiv] = true;
98:
99: #endregion
100:
101: }
102:
103: public DLASY2()
104: {
105:
106:
107: #region Dependencies (Initialization)
108:
109: IDAMAX idamax = new IDAMAX();
110: LSAME lsame = new LSAME();
111: DLAMC3 dlamc3 = new DLAMC3();
112: DCOPY dcopy = new DCOPY();
113: DSWAP dswap = new DSWAP();
114: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
115: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
116: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
117: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
118: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
119:
120: #endregion
121:
122:
123: #region Set Dependencies
124:
125: this._idamax = idamax; this._dlamch = dlamch; this._dcopy = dcopy; this._dswap = dswap;
126:
127: #endregion
128:
129:
130: #region Data Inicializacion
131:
132: //LOCU12/3,4,1,2
133: LOCU12[1 + o_locu12] = 3;
134: LOCU12[2 + o_locu12] = 4;
135: LOCU12[3 + o_locu12] = 1;
136: LOCU12[4 + o_locu12] = 2;
137: //LOCL21/2,1,4,3
138: LOCL21[1 + o_locl21] = 2;
139: LOCL21[2 + o_locl21] = 1;
140: LOCL21[3 + o_locl21] = 4;
141: LOCL21[4 + o_locl21] = 3;
142: //LOCU22/4,3,2,1
143: LOCU22[1 + o_locu22] = 4;
144: LOCU22[2 + o_locu22] = 3;
145: LOCU22[3 + o_locu22] = 2;
146: LOCU22[4 + o_locu22] = 1;
147: //XSWPIV/.FALSE.,.FALSE.,.TRUE.,.TRUE.
148: XSWPIV[1 + o_xswpiv] = false;
149: XSWPIV[2 + o_xswpiv] = false;
150: XSWPIV[3 + o_xswpiv] = true;
151: XSWPIV[4 + o_xswpiv] = true;
152: //BSWPIV/.FALSE.,.TRUE.,.FALSE.,.TRUE.
153: BSWPIV[1 + o_bswpiv] = false;
154: BSWPIV[2 + o_bswpiv] = true;
155: BSWPIV[3 + o_bswpiv] = false;
156: BSWPIV[4 + o_bswpiv] = true;
157:
158: #endregion
159:
160: }
161: /// <summary>
162: /// Purpose
163: /// =======
164: ///
165: /// DLASY2 solves for the N1 by N2 matrix X, 1 .LE. N1,N2 .LE. 2, in
166: ///
167: /// op(TL)*X + ISGN*X*op(TR) = SCALE*B,
168: ///
169: /// where TL is N1 by N1, TR is N2 by N2, B is N1 by N2, and ISGN = 1 or
170: /// -1. op(T) = T or T', where T' denotes the transpose of T.
171: ///
172: ///</summary>
173: /// <param name="LTRANL">
174: /// (input) LOGICAL
175: /// On entry, LTRANL specifies the op(TL):
176: /// = .FALSE., op(TL) = TL,
177: /// = .TRUE., op(TL) = TL'.
178: ///</param>
179: /// <param name="LTRANR">
180: /// (input) LOGICAL
181: /// On entry, LTRANR specifies the op(TR):
182: /// = .FALSE., op(TR) = TR,
183: /// = .TRUE., op(TR) = TR'.
184: ///</param>
185: /// <param name="ISGN">
186: /// (input) INTEGER
187: /// On entry, ISGN specifies the sign of the equation
188: /// as described before. ISGN may only be 1 or -1.
189: ///</param>
190: /// <param name="N1">
191: /// (input) INTEGER
192: /// On entry, N1 specifies the order of matrix TL.
193: /// N1 may only be 0, 1 or 2.
194: ///</param>
195: /// <param name="N2">
196: /// (input) INTEGER
197: /// On entry, N2 specifies the order of matrix TR.
198: /// N2 may only be 0, 1 or 2.
199: ///</param>
200: /// <param name="TL">
201: /// (input) DOUBLE PRECISION array, dimension (LDTL,2)
202: /// On entry, TL contains an N1 by N1 matrix.
203: ///</param>
204: /// <param name="LDTL">
205: /// (input) INTEGER
206: /// The leading dimension of the matrix TL. LDTL .GE. max(1,N1).
207: ///</param>
208: /// <param name="TR">
209: /// (input) DOUBLE PRECISION array, dimension (LDTR,2)
210: /// On entry, TR contains an N2 by N2 matrix.
211: ///</param>
212: /// <param name="LDTR">
213: /// (input) INTEGER
214: /// The leading dimension of the matrix TR. LDTR .GE. max(1,N2).
215: ///</param>
216: /// <param name="B">
217: /// (input) DOUBLE PRECISION array, dimension (LDB,2)
218: /// On entry, the N1 by N2 matrix B contains the right-hand
219: /// side of the equation.
220: ///</param>
221: /// <param name="LDB">
222: /// (input) INTEGER
223: /// The leading dimension of the matrix B. LDB .GE. max(1,N1).
224: ///</param>
225: /// <param name="SCALE">
226: /// (output) DOUBLE PRECISION
227: /// On exit, SCALE contains the scale factor. SCALE is chosen
228: /// less than or equal to 1 to prevent the solution overflowing.
229: ///</param>
230: /// <param name="X">
231: /// (output) DOUBLE PRECISION array, dimension (LDX,2)
232: /// On exit, X contains the N1 by N2 solution.
233: ///</param>
234: /// <param name="LDX">
235: /// (input) INTEGER
236: /// The leading dimension of the matrix X. LDX .GE. max(1,N1).
237: ///</param>
238: /// <param name="XNORM">
239: /// (output) DOUBLE PRECISION
240: /// On exit, XNORM is the infinity-norm of the solution.
241: ///</param>
242: /// <param name="INFO">
243: /// (output) INTEGER
244: /// On exit, INFO is set to
245: /// 0: successful exit.
246: /// 1: TL and TR have too close eigenvalues, so TL or
247: /// TR is perturbed to get a nonsingular equation.
248: /// NOTE: In the interests of speed, this routine does not
249: /// check the inputs for errors.
250: ///</param>
251: public void Run(bool LTRANL, bool LTRANR, int ISGN, int N1, int N2, double[] TL, int offset_tl
252: , int LDTL, double[] TR, int offset_tr, int LDTR, double[] B, int offset_b, int LDB, ref double SCALE
253: , ref double[] X, int offset_x, int LDX, ref double XNORM, ref int INFO)
254: {
255:
256: #region Array Index Correction
257:
258: int o_tl = -1 - LDTL + offset_tl; int o_tr = -1 - LDTR + offset_tr; int o_b = -1 - LDB + offset_b;
259: int o_x = -1 - LDX + offset_x;
260:
261: #endregion
262:
263:
264: #region Prolog
265:
266: // *
267: // * -- LAPACK auxiliary routine (version 3.1) --
268: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
269: // * November 2006
270: // *
271: // * .. Scalar Arguments ..
272: // * ..
273: // * .. Array Arguments ..
274: // * ..
275: // *
276: // * Purpose
277: // * =======
278: // *
279: // * DLASY2 solves for the N1 by N2 matrix X, 1 <= N1,N2 <= 2, in
280: // *
281: // * op(TL)*X + ISGN*X*op(TR) = SCALE*B,
282: // *
283: // * where TL is N1 by N1, TR is N2 by N2, B is N1 by N2, and ISGN = 1 or
284: // * -1. op(T) = T or T', where T' denotes the transpose of T.
285: // *
286: // * Arguments
287: // * =========
288: // *
289: // * LTRANL (input) LOGICAL
290: // * On entry, LTRANL specifies the op(TL):
291: // * = .FALSE., op(TL) = TL,
292: // * = .TRUE., op(TL) = TL'.
293: // *
294: // * LTRANR (input) LOGICAL
295: // * On entry, LTRANR specifies the op(TR):
296: // * = .FALSE., op(TR) = TR,
297: // * = .TRUE., op(TR) = TR'.
298: // *
299: // * ISGN (input) INTEGER
300: // * On entry, ISGN specifies the sign of the equation
301: // * as described before. ISGN may only be 1 or -1.
302: // *
303: // * N1 (input) INTEGER
304: // * On entry, N1 specifies the order of matrix TL.
305: // * N1 may only be 0, 1 or 2.
306: // *
307: // * N2 (input) INTEGER
308: // * On entry, N2 specifies the order of matrix TR.
309: // * N2 may only be 0, 1 or 2.
310: // *
311: // * TL (input) DOUBLE PRECISION array, dimension (LDTL,2)
312: // * On entry, TL contains an N1 by N1 matrix.
313: // *
314: // * LDTL (input) INTEGER
315: // * The leading dimension of the matrix TL. LDTL >= max(1,N1).
316: // *
317: // * TR (input) DOUBLE PRECISION array, dimension (LDTR,2)
318: // * On entry, TR contains an N2 by N2 matrix.
319: // *
320: // * LDTR (input) INTEGER
321: // * The leading dimension of the matrix TR. LDTR >= max(1,N2).
322: // *
323: // * B (input) DOUBLE PRECISION array, dimension (LDB,2)
324: // * On entry, the N1 by N2 matrix B contains the right-hand
325: // * side of the equation.
326: // *
327: // * LDB (input) INTEGER
328: // * The leading dimension of the matrix B. LDB >= max(1,N1).
329: // *
330: // * SCALE (output) DOUBLE PRECISION
331: // * On exit, SCALE contains the scale factor. SCALE is chosen
332: // * less than or equal to 1 to prevent the solution overflowing.
333: // *
334: // * X (output) DOUBLE PRECISION array, dimension (LDX,2)
335: // * On exit, X contains the N1 by N2 solution.
336: // *
337: // * LDX (input) INTEGER
338: // * The leading dimension of the matrix X. LDX >= max(1,N1).
339: // *
340: // * XNORM (output) DOUBLE PRECISION
341: // * On exit, XNORM is the infinity-norm of the solution.
342: // *
343: // * INFO (output) INTEGER
344: // * On exit, INFO is set to
345: // * 0: successful exit.
346: // * 1: TL and TR have too close eigenvalues, so TL or
347: // * TR is perturbed to get a nonsingular equation.
348: // * NOTE: In the interests of speed, this routine does not
349: // * check the inputs for errors.
350: // *
351: // * =====================================================================
352: // *
353: // * .. Parameters ..
354: // * ..
355: // * .. Local Scalars ..
356: // * ..
357: // * .. Local Arrays ..
358: // * ..
359: // * .. External Functions ..
360: // * ..
361: // * .. External Subroutines ..
362: // * ..
363: // * .. Intrinsic Functions ..
364: // INTRINSIC ABS, MAX;
365: // * ..
366: // * .. Data statements ..
367: // * ..
368: // * .. Executable Statements ..
369: // *
370: // * Do not check the input parameters for errors
371: // *
372:
373: #endregion
374:
375:
376: #region Body
377:
378: INFO = 0;
379: // *
380: // * Quick return if possible
381: // *
382: if (N1 == 0 || N2 == 0) return;
383: // *
384: // * Set constants to control overflow
385: // *
386: EPS = this._dlamch.Run("P");
387: SMLNUM = this._dlamch.Run("S") / EPS;
388: SGN = ISGN;
389: // *
390: K = N1 + N1 + N2 - 2;
391: switch (K)
392: {
393: case 1: goto LABEL10;
394: case 2: goto LABEL20;
395: case 3: goto LABEL30;
396: case 4: goto LABEL50;
397: }
398: // *
399: // * 1 by 1: TL11*X + SGN*X*TR11 = B11
400: // *
401: LABEL10:;
402: TAU1 = TL[1+1 * LDTL + o_tl] + SGN * TR[1+1 * LDTR + o_tr];
403: BET = Math.Abs(TAU1);
404: if (BET <= SMLNUM)
405: {
406: TAU1 = SMLNUM;
407: BET = SMLNUM;
408: INFO = 1;
409: }
410: // *
411: SCALE = ONE;
412: GAM = Math.Abs(B[1+1 * LDB + o_b]);
413: if (SMLNUM * GAM > BET) SCALE = ONE / GAM;
414: // *
415: X[1+1 * LDX + o_x] = (B[1+1 * LDB + o_b] * SCALE) / TAU1;
416: XNORM = Math.Abs(X[1+1 * LDX + o_x]);
417: return;
418: // *
419: // * 1 by 2:
420: // * TL11*[X11 X12] + ISGN*[X11 X12]*op[TR11 TR12] = [B11 B12]
421: // * [TR21 TR22]
422: // *
423: LABEL20:;
424: // *
425: SMIN = Math.Max(EPS * Math.Max(Math.Abs(TL[1+1 * LDTL + o_tl]), Math.Max(Math.Abs(TR[1+1 * LDTR + o_tr]), Math.Max(Math.Abs(TR[1+2 * LDTR + o_tr]), Math.Max(Math.Abs(TR[2+1 * LDTR + o_tr]), Math.Abs(TR[2+2 * LDTR + o_tr]))))), SMLNUM);
426: TMP[1 + o_tmp] = TL[1+1 * LDTL + o_tl] + SGN * TR[1+1 * LDTR + o_tr];
427: TMP[4 + o_tmp] = TL[1+1 * LDTL + o_tl] + SGN * TR[2+2 * LDTR + o_tr];
428: if (LTRANR)
429: {
430: TMP[2 + o_tmp] = SGN * TR[2+1 * LDTR + o_tr];
431: TMP[3 + o_tmp] = SGN * TR[1+2 * LDTR + o_tr];
432: }
433: else
434: {
435: TMP[2 + o_tmp] = SGN * TR[1+2 * LDTR + o_tr];
436: TMP[3 + o_tmp] = SGN * TR[2+1 * LDTR + o_tr];
437: }
438: BTMP[1 + o_btmp] = B[1+1 * LDB + o_b];
439: BTMP[2 + o_btmp] = B[1+2 * LDB + o_b];
440: goto LABEL40;
441: // *
442: // * 2 by 1:
443: // * op[TL11 TL12]*[X11] + ISGN* [X11]*TR11 = [B11]
444: // * [TL21 TL22] [X21] [X21] [B21]
445: // *
446: LABEL30:;
447: SMIN = Math.Max(EPS * Math.Max(Math.Abs(TR[1+1 * LDTR + o_tr]), Math.Max(Math.Abs(TL[1+1 * LDTL + o_tl]), Math.Max(Math.Abs(TL[1+2 * LDTL + o_tl]), Math.Max(Math.Abs(TL[2+1 * LDTL + o_tl]), Math.Abs(TL[2+2 * LDTL + o_tl]))))), SMLNUM);
448: TMP[1 + o_tmp] = TL[1+1 * LDTL + o_tl] + SGN * TR[1+1 * LDTR + o_tr];
449: TMP[4 + o_tmp] = TL[2+2 * LDTL + o_tl] + SGN * TR[1+1 * LDTR + o_tr];
450: if (LTRANL)
451: {
452: TMP[2 + o_tmp] = TL[1+2 * LDTL + o_tl];
453: TMP[3 + o_tmp] = TL[2+1 * LDTL + o_tl];
454: }
455: else
456: {
457: TMP[2 + o_tmp] = TL[2+1 * LDTL + o_tl];
458: TMP[3 + o_tmp] = TL[1+2 * LDTL + o_tl];
459: }
460: BTMP[1 + o_btmp] = B[1+1 * LDB + o_b];
461: BTMP[2 + o_btmp] = B[2+1 * LDB + o_b];
462: LABEL40:;
463: // *
464: // * Solve 2 by 2 system using complete pivoting.
465: // * Set pivots less than SMIN to SMIN.
466: // *
467: IPIV = this._idamax.Run(4, TMP, offset_tmp, 1);
468: U11 = TMP[IPIV + o_tmp];
469: if (Math.Abs(U11) <= SMIN)
470: {
471: INFO = 1;
472: U11 = SMIN;
473: }
474: U12 = TMP[LOCU12[IPIV + o_locu12] + o_tmp];
475: L21 = TMP[LOCL21[IPIV + o_locl21] + o_tmp] / U11;
476: U22 = TMP[LOCU22[IPIV + o_locu22] + o_tmp] - U12 * L21;
477: XSWAP = XSWPIV[IPIV + o_xswpiv];
478: BSWAP = BSWPIV[IPIV + o_bswpiv];
479: if (Math.Abs(U22) <= SMIN)
480: {
481: INFO = 1;
482: U22 = SMIN;
483: }
484: if (BSWAP)
485: {
486: TEMP = BTMP[2 + o_btmp];
487: BTMP[2 + o_btmp] = BTMP[1 + o_btmp] - L21 * TEMP;
488: BTMP[1 + o_btmp] = TEMP;
489: }
490: else
491: {
492: BTMP[2 + o_btmp] = BTMP[2 + o_btmp] - L21 * BTMP[1 + o_btmp];
493: }
494: SCALE = ONE;
495: if ((TWO * SMLNUM) * Math.Abs(BTMP[2 + o_btmp]) > Math.Abs(U22) || (TWO * SMLNUM) * Math.Abs(BTMP[1 + o_btmp]) > Math.Abs(U11))
496: {
497: SCALE = HALF / Math.Max(Math.Abs(BTMP[1 + o_btmp]), Math.Abs(BTMP[2 + o_btmp]));
498: BTMP[1 + o_btmp] = BTMP[1 + o_btmp] * SCALE;
499: BTMP[2 + o_btmp] = BTMP[2 + o_btmp] * SCALE;
500: }
501: X2[2 + o_x2] = BTMP[2 + o_btmp] / U22;
502: X2[1 + o_x2] = BTMP[1 + o_btmp] / U11 - (U12 / U11) * X2[2 + o_x2];
503: if (XSWAP)
504: {
505: TEMP = X2[2 + o_x2];
506: X2[2 + o_x2] = X2[1 + o_x2];
507: X2[1 + o_x2] = TEMP;
508: }
509: X[1+1 * LDX + o_x] = X2[1 + o_x2];
510: if (N1 == 1)
511: {
512: X[1+2 * LDX + o_x] = X2[2 + o_x2];
513: XNORM = Math.Abs(X[1+1 * LDX + o_x]) + Math.Abs(X[1+2 * LDX + o_x]);
514: }
515: else
516: {
517: X[2+1 * LDX + o_x] = X2[2 + o_x2];
518: XNORM = Math.Max(Math.Abs(X[1+1 * LDX + o_x]), Math.Abs(X[2+1 * LDX + o_x]));
519: }
520: return;
521: // *
522: // * 2 by 2:
523: // * op[TL11 TL12]*[X11 X12] +ISGN* [X11 X12]*op[TR11 TR12] = [B11 B12]
524: // * [TL21 TL22] [X21 X22] [X21 X22] [TR21 TR22] [B21 B22]
525: // *
526: // * Solve equivalent 4 by 4 system using complete pivoting.
527: // * Set pivots less than SMIN to SMIN.
528: // *
529: LABEL50:;
530: SMIN = Math.Max(Math.Abs(TR[1+1 * LDTR + o_tr]), Math.Max(Math.Abs(TR[1+2 * LDTR + o_tr]), Math.Max(Math.Abs(TR[2+1 * LDTR + o_tr]), Math.Abs(TR[2+2 * LDTR + o_tr]))));
531: SMIN = Math.Max(SMIN, Math.Max(Math.Abs(TL[1+1 * LDTL + o_tl]), Math.Max(Math.Abs(TL[1+2 * LDTL + o_tl]), Math.Max(Math.Abs(TL[2+1 * LDTL + o_tl]), Math.Abs(TL[2+2 * LDTL + o_tl])))));
532: SMIN = Math.Max(EPS * SMIN, SMLNUM);
533: BTMP[1 + o_btmp] = ZERO;
534: this._dcopy.Run(16, BTMP, offset_btmp, 0, ref T16, offset_t16, 1);
535: T16[1+1 * 4 + o_t16] = TL[1+1 * LDTL + o_tl] + SGN * TR[1+1 * LDTR + o_tr];
536: T16[2+2 * 4 + o_t16] = TL[2+2 * LDTL + o_tl] + SGN * TR[1+1 * LDTR + o_tr];
537: T16[3+3 * 4 + o_t16] = TL[1+1 * LDTL + o_tl] + SGN * TR[2+2 * LDTR + o_tr];
538: T16[4+4 * 4 + o_t16] = TL[2+2 * LDTL + o_tl] + SGN * TR[2+2 * LDTR + o_tr];
539: if (LTRANL)
540: {
541: T16[1+2 * 4 + o_t16] = TL[2+1 * LDTL + o_tl];
542: T16[2+1 * 4 + o_t16] = TL[1+2 * LDTL + o_tl];
543: T16[3+4 * 4 + o_t16] = TL[2+1 * LDTL + o_tl];
544: T16[4+3 * 4 + o_t16] = TL[1+2 * LDTL + o_tl];
545: }
546: else
547: {
548: T16[1+2 * 4 + o_t16] = TL[1+2 * LDTL + o_tl];
549: T16[2+1 * 4 + o_t16] = TL[2+1 * LDTL + o_tl];
550: T16[3+4 * 4 + o_t16] = TL[1+2 * LDTL + o_tl];
551: T16[4+3 * 4 + o_t16] = TL[2+1 * LDTL + o_tl];
552: }
553: if (LTRANR)
554: {
555: T16[1+3 * 4 + o_t16] = SGN * TR[1+2 * LDTR + o_tr];
556: T16[2+4 * 4 + o_t16] = SGN * TR[1+2 * LDTR + o_tr];
557: T16[3+1 * 4 + o_t16] = SGN * TR[2+1 * LDTR + o_tr];
558: T16[4+2 * 4 + o_t16] = SGN * TR[2+1 * LDTR + o_tr];
559: }
560: else
561: {
562: T16[1+3 * 4 + o_t16] = SGN * TR[2+1 * LDTR + o_tr];
563: T16[2+4 * 4 + o_t16] = SGN * TR[2+1 * LDTR + o_tr];
564: T16[3+1 * 4 + o_t16] = SGN * TR[1+2 * LDTR + o_tr];
565: T16[4+2 * 4 + o_t16] = SGN * TR[1+2 * LDTR + o_tr];
566: }
567: BTMP[1 + o_btmp] = B[1+1 * LDB + o_b];
568: BTMP[2 + o_btmp] = B[2+1 * LDB + o_b];
569: BTMP[3 + o_btmp] = B[1+2 * LDB + o_b];
570: BTMP[4 + o_btmp] = B[2+2 * LDB + o_b];
571: // *
572: // * Perform elimination
573: // *
574: for (I = 1; I <= 3; I++)
575: {
576: XMAX = ZERO;
577: for (IP = I; IP <= 4; IP++)
578: {
579: for (JP = I; JP <= 4; JP++)
580: {
581: if (Math.Abs(T16[IP+JP * 4 + o_t16]) >= XMAX)
582: {
583: XMAX = Math.Abs(T16[IP+JP * 4 + o_t16]);
584: IPSV = IP;
585: JPSV = JP;
586: }
587: }
588: }
589: if (IPSV != I)
590: {
591: this._dswap.Run(4, ref T16, IPSV+1 * 4 + o_t16, 4, ref T16, I+1 * 4 + o_t16, 4);
592: TEMP = BTMP[I + o_btmp];
593: BTMP[I + o_btmp] = BTMP[IPSV + o_btmp];
594: BTMP[IPSV + o_btmp] = TEMP;
595: }
596: if (JPSV != I) this._dswap.Run(4, ref T16, 1+JPSV * 4 + o_t16, 1, ref T16, 1+I * 4 + o_t16, 1);
597: JPIV[I + o_jpiv] = JPSV;
598: if (Math.Abs(T16[I+I * 4 + o_t16]) < SMIN)
599: {
600: INFO = 1;
601: T16[I+I * 4 + o_t16] = SMIN;
602: }
603: for (J = I + 1; J <= 4; J++)
604: {
605: T16[J+I * 4 + o_t16] = T16[J+I * 4 + o_t16] / T16[I+I * 4 + o_t16];
606: BTMP[J + o_btmp] = BTMP[J + o_btmp] - T16[J+I * 4 + o_t16] * BTMP[I + o_btmp];
607: for (K = I + 1; K <= 4; K++)
608: {
609: T16[J+K * 4 + o_t16] = T16[J+K * 4 + o_t16] - T16[J+I * 4 + o_t16] * T16[I+K * 4 + o_t16];
610: }
611: }
612: }
613: if (Math.Abs(T16[4+4 * 4 + o_t16]) < SMIN) T16[4+4 * 4 + o_t16] = SMIN;
614: SCALE = ONE;
615: if ((EIGHT * SMLNUM) * Math.Abs(BTMP[1 + o_btmp]) > Math.Abs(T16[1+1 * 4 + o_t16]) || (EIGHT * SMLNUM) * Math.Abs(BTMP[2 + o_btmp]) > Math.Abs(T16[2+2 * 4 + o_t16]) || (EIGHT * SMLNUM) * Math.Abs(BTMP[3 + o_btmp]) > Math.Abs(T16[3+3 * 4 + o_t16]) || (EIGHT * SMLNUM) * Math.Abs(BTMP[4 + o_btmp]) > Math.Abs(T16[4+4 * 4 + o_t16]))
616: {
617: SCALE = (ONE / EIGHT) / Math.Max(Math.Abs(BTMP[1 + o_btmp]), Math.Max(Math.Abs(BTMP[2 + o_btmp]), Math.Max(Math.Abs(BTMP[3 + o_btmp]), Math.Abs(BTMP[4 + o_btmp]))));
618: BTMP[1 + o_btmp] = BTMP[1 + o_btmp] * SCALE;
619: BTMP[2 + o_btmp] = BTMP[2 + o_btmp] * SCALE;
620: BTMP[3 + o_btmp] = BTMP[3 + o_btmp] * SCALE;
621: BTMP[4 + o_btmp] = BTMP[4 + o_btmp] * SCALE;
622: }
623: for (I = 1; I <= 4; I++)
624: {
625: K = 5 - I;
626: TEMP = ONE / T16[K+K * 4 + o_t16];
627: TMP[K + o_tmp] = BTMP[K + o_btmp] * TEMP;
628: for (J = K + 1; J <= 4; J++)
629: {
630: TMP[K + o_tmp] = TMP[K + o_tmp] - (TEMP * T16[K+J * 4 + o_t16]) * TMP[J + o_tmp];
631: }
632: }
633: for (I = 1; I <= 3; I++)
634: {
635: if (JPIV[4 - I + o_jpiv] != 4 - I)
636: {
637: TEMP = TMP[4 - I + o_tmp];
638: TMP[4 - I + o_tmp] = TMP[JPIV[4 - I + o_jpiv] + o_tmp];
639: TMP[JPIV[4 - I + o_jpiv] + o_tmp] = TEMP;
640: }
641: }
642: X[1+1 * LDX + o_x] = TMP[1 + o_tmp];
643: X[2+1 * LDX + o_x] = TMP[2 + o_tmp];
644: X[1+2 * LDX + o_x] = TMP[3 + o_tmp];
645: X[2+2 * LDX + o_x] = TMP[4 + o_tmp];
646: XNORM = Math.Max(Math.Abs(TMP[1 + o_tmp]) + Math.Abs(TMP[3 + o_tmp]), Math.Abs(TMP[2 + o_tmp]) + Math.Abs(TMP[4 + o_tmp]));
647: return;
648: // *
649: // * End of DLASY2
650: // *
651:
652: #endregion
653:
654: }
655: }
656: }