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.0) --
21: /// Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
22: /// Courant Institute, Argonne National Lab, and Rice University
23: /// June 30, 1992
24: /// Purpose
25: /// =======
26: ///
27: /// DLATRS solves one of the triangular systems
28: ///
29: /// A *x = s*b or A'*x = s*b
30: ///
31: /// with scaling to prevent overflow. Here A is an upper or lower
32: /// triangular matrix, A' denotes the transpose of A, x and b are
33: /// n-element vectors, and s is a scaling factor, usually less than
34: /// or equal to 1, chosen so that the components of x will be less than
35: /// the overflow threshold. If the unscaled problem will not cause
36: /// overflow, the Level 2 BLAS routine DTRSV is called. If the matrix A
37: /// is singular (A(j,j) = 0 for some j), then s is set to 0 and a
38: /// non-trivial solution to A*x = 0 is returned.
39: ///
40: ///</summary>
41: public class DLATRS
42: {
43:
44:
45: #region Dependencies
46:
47: LSAME _lsame; IDAMAX _idamax; DASUM _dasum; DDOT _ddot; DLAMCH _dlamch; DAXPY _daxpy; DSCAL _dscal; DTRSV _dtrsv;
48: XERBLA _xerbla;
49:
50: #endregion
51:
52:
53: #region Fields
54:
55: const double ZERO = 0.0E+0; const double HALF = 0.5E+0; const double ONE = 1.0E+0; bool NOTRAN = false;
56: bool NOUNIT = false;bool UPPER = false; int I = 0; int IMAX = 0; int J = 0; int JFIRST = 0; int JINC = 0; int JLAST = 0;
57: double BIGNUM = 0;double GROW = 0; double REC = 0; double SMLNUM = 0; double SUMJ = 0; double TJJ = 0; double TJJS = 0;
58: double TMAX = 0;double TSCAL = 0; double USCAL = 0; double XBND = 0; double XJ = 0; double XMAX = 0;
59:
60: #endregion
61:
62: public DLATRS(LSAME lsame, IDAMAX idamax, DASUM dasum, DDOT ddot, DLAMCH dlamch, DAXPY daxpy, DSCAL dscal, DTRSV dtrsv, XERBLA xerbla)
63: {
64:
65:
66: #region Set Dependencies
67:
68: this._lsame = lsame; this._idamax = idamax; this._dasum = dasum; this._ddot = ddot; this._dlamch = dlamch;
69: this._daxpy = daxpy;this._dscal = dscal; this._dtrsv = dtrsv; this._xerbla = xerbla;
70:
71: #endregion
72:
73: }
74:
75: public DLATRS()
76: {
77:
78:
79: #region Dependencies (Initialization)
80:
81: LSAME lsame = new LSAME();
82: IDAMAX idamax = new IDAMAX();
83: DASUM dasum = new DASUM();
84: DDOT ddot = new DDOT();
85: DLAMC3 dlamc3 = new DLAMC3();
86: DAXPY daxpy = new DAXPY();
87: DSCAL dscal = new DSCAL();
88: XERBLA xerbla = new XERBLA();
89: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
90: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
91: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
92: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
93: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
94: DTRSV dtrsv = new DTRSV(lsame, xerbla);
95:
96: #endregion
97:
98:
99: #region Set Dependencies
100:
101: this._lsame = lsame; this._idamax = idamax; this._dasum = dasum; this._ddot = ddot; this._dlamch = dlamch;
102: this._daxpy = daxpy;this._dscal = dscal; this._dtrsv = dtrsv; this._xerbla = xerbla;
103:
104: #endregion
105:
106: }
107: /// <summary>
108: /// Purpose
109: /// =======
110: ///
111: /// DLATRS solves one of the triangular systems
112: ///
113: /// A *x = s*b or A'*x = s*b
114: ///
115: /// with scaling to prevent overflow. Here A is an upper or lower
116: /// triangular matrix, A' denotes the transpose of A, x and b are
117: /// n-element vectors, and s is a scaling factor, usually less than
118: /// or equal to 1, chosen so that the components of x will be less than
119: /// the overflow threshold. If the unscaled problem will not cause
120: /// overflow, the Level 2 BLAS routine DTRSV is called. If the matrix A
121: /// is singular (A(j,j) = 0 for some j), then s is set to 0 and a
122: /// non-trivial solution to A*x = 0 is returned.
123: ///
124: ///</summary>
125: /// <param name="UPLO">
126: /// (input) CHARACTER*1
127: /// Specifies whether the matrix A is upper or lower triangular.
128: /// = 'U': Upper triangular
129: /// = 'L': Lower triangular
130: ///</param>
131: /// <param name="TRANS">
132: /// (input) CHARACTER*1
133: /// Specifies the operation applied to A.
134: /// = 'N': Solve A * x = s*b (No transpose)
135: /// = 'T': Solve A'* x = s*b (Transpose)
136: /// = 'C': Solve A'* x = s*b (Conjugate transpose = Transpose)
137: ///</param>
138: /// <param name="DIAG">
139: /// (input) CHARACTER*1
140: /// Specifies whether or not the matrix A is unit triangular.
141: /// = 'N': Non-unit triangular
142: /// = 'U': Unit triangular
143: ///</param>
144: /// <param name="NORMIN">
145: /// (input) CHARACTER*1
146: /// Specifies whether CNORM has been set or not.
147: /// = 'Y': CNORM contains the column norms on entry
148: /// = 'N': CNORM is not set on entry. On exit, the norms will
149: /// be computed and stored in CNORM.
150: ///</param>
151: /// <param name="N">
152: /// (input) INTEGER
153: /// The order of the matrix A. N .GE. 0.
154: ///</param>
155: /// <param name="A">
156: /// *x = s*b or A'*x = s*b
157: ///</param>
158: /// <param name="LDA">
159: /// (input) INTEGER
160: /// The leading dimension of the array A. LDA .GE. max (1,N).
161: ///</param>
162: /// <param name="X">
163: /// (input/output) DOUBLE PRECISION array, dimension (N)
164: /// On entry, the right hand side b of the triangular system.
165: /// On exit, X is overwritten by the solution vector x.
166: ///</param>
167: /// <param name="SCALE">
168: /// (output) DOUBLE PRECISION
169: /// The scaling factor s for the triangular system
170: /// A * x = s*b or A'* x = s*b.
171: /// If SCALE = 0, the matrix A is singular or badly scaled, and
172: /// the vector x is an exact or approximate solution to A*x = 0.
173: ///</param>
174: /// <param name="CNORM">
175: /// (input or output) DOUBLE PRECISION array, dimension (N)
176: ///
177: /// If NORMIN = 'Y', CNORM is an input argument and CNORM(j)
178: /// contains the norm of the off-diagonal part of the j-th column
179: /// of A. If TRANS = 'N', CNORM(j) must be greater than or equal
180: /// to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)
181: /// must be greater than or equal to the 1-norm.
182: ///
183: /// If NORMIN = 'N', CNORM is an output argument and CNORM(j)
184: /// returns the 1-norm of the offdiagonal part of the j-th column
185: /// of A.
186: ///</param>
187: /// <param name="INFO">
188: /// (output) INTEGER
189: /// = 0: successful exit
190: /// .LT. 0: if INFO = -k, the k-th argument had an illegal value
191: ///</param>
192: public void Run(string UPLO, string TRANS, string DIAG, string NORMIN, int N, double[] A, int offset_a
193: , int LDA, ref double[] X, int offset_x, ref double SCALE, ref double[] CNORM, int offset_cnorm, ref int INFO)
194: {
195:
196: #region Array Index Correction
197:
198: int o_a = -1 - LDA + offset_a; int o_x = -1 + offset_x; int o_cnorm = -1 + offset_cnorm;
199:
200: #endregion
201:
202:
203: #region Strings
204:
205: UPLO = UPLO.Substring(0, 1); TRANS = TRANS.Substring(0, 1); DIAG = DIAG.Substring(0, 1);
206: NORMIN = NORMIN.Substring(0, 1);
207:
208: #endregion
209:
210:
211: #region Prolog
212:
213: // *
214: // * -- LAPACK auxiliary routine (version 3.0) --
215: // * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
216: // * Courant Institute, Argonne National Lab, and Rice University
217: // * June 30, 1992
218: // *
219: // * .. Scalar Arguments ..
220: // * ..
221: // * .. Array Arguments ..
222: // * ..
223: // *
224: // * Purpose
225: // * =======
226: // *
227: // * DLATRS solves one of the triangular systems
228: // *
229: // * A *x = s*b or A'*x = s*b
230: // *
231: // * with scaling to prevent overflow. Here A is an upper or lower
232: // * triangular matrix, A' denotes the transpose of A, x and b are
233: // * n-element vectors, and s is a scaling factor, usually less than
234: // * or equal to 1, chosen so that the components of x will be less than
235: // * the overflow threshold. If the unscaled problem will not cause
236: // * overflow, the Level 2 BLAS routine DTRSV is called. If the matrix A
237: // * is singular (A(j,j) = 0 for some j), then s is set to 0 and a
238: // * non-trivial solution to A*x = 0 is returned.
239: // *
240: // * Arguments
241: // * =========
242: // *
243: // * UPLO (input) CHARACTER*1
244: // * Specifies whether the matrix A is upper or lower triangular.
245: // * = 'U': Upper triangular
246: // * = 'L': Lower triangular
247: // *
248: // * TRANS (input) CHARACTER*1
249: // * Specifies the operation applied to A.
250: // * = 'N': Solve A * x = s*b (No transpose)
251: // * = 'T': Solve A'* x = s*b (Transpose)
252: // * = 'C': Solve A'* x = s*b (Conjugate transpose = Transpose)
253: // *
254: // * DIAG (input) CHARACTER*1
255: // * Specifies whether or not the matrix A is unit triangular.
256: // * = 'N': Non-unit triangular
257: // * = 'U': Unit triangular
258: // *
259: // * NORMIN (input) CHARACTER*1
260: // * Specifies whether CNORM has been set or not.
261: // * = 'Y': CNORM contains the column norms on entry
262: // * = 'N': CNORM is not set on entry. On exit, the norms will
263: // * be computed and stored in CNORM.
264: // *
265: // * N (input) INTEGER
266: // * The order of the matrix A. N >= 0.
267: // *
268: // * A (input) DOUBLE PRECISION array, dimension (LDA,N)
269: // * The triangular matrix A. If UPLO = 'U', the leading n by n
270: // * upper triangular part of the array A contains the upper
271: // * triangular matrix, and the strictly lower triangular part of
272: // * A is not referenced. If UPLO = 'L', the leading n by n lower
273: // * triangular part of the array A contains the lower triangular
274: // * matrix, and the strictly upper triangular part of A is not
275: // * referenced. If DIAG = 'U', the diagonal elements of A are
276: // * also not referenced and are assumed to be 1.
277: // *
278: // * LDA (input) INTEGER
279: // * The leading dimension of the array A. LDA >= max (1,N).
280: // *
281: // * X (input/output) DOUBLE PRECISION array, dimension (N)
282: // * On entry, the right hand side b of the triangular system.
283: // * On exit, X is overwritten by the solution vector x.
284: // *
285: // * SCALE (output) DOUBLE PRECISION
286: // * The scaling factor s for the triangular system
287: // * A * x = s*b or A'* x = s*b.
288: // * If SCALE = 0, the matrix A is singular or badly scaled, and
289: // * the vector x is an exact or approximate solution to A*x = 0.
290: // *
291: // * CNORM (input or output) DOUBLE PRECISION array, dimension (N)
292: // *
293: // * If NORMIN = 'Y', CNORM is an input argument and CNORM(j)
294: // * contains the norm of the off-diagonal part of the j-th column
295: // * of A. If TRANS = 'N', CNORM(j) must be greater than or equal
296: // * to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)
297: // * must be greater than or equal to the 1-norm.
298: // *
299: // * If NORMIN = 'N', CNORM is an output argument and CNORM(j)
300: // * returns the 1-norm of the offdiagonal part of the j-th column
301: // * of A.
302: // *
303: // * INFO (output) INTEGER
304: // * = 0: successful exit
305: // * < 0: if INFO = -k, the k-th argument had an illegal value
306: // *
307: // * Further Details
308: // * ======= =======
309: // *
310: // * A rough bound on x is computed; if that is less than overflow, DTRSV
311: // * is called, otherwise, specific code is used which checks for possible
312: // * overflow or divide-by-zero at every operation.
313: // *
314: // * A columnwise scheme is used for solving A*x = b. The basic algorithm
315: // * if A is lower triangular is
316: // *
317: // * x[1:n] := b[1:n]
318: // * for j = 1, ..., n
319: // * x(j) := x(j) / A(j,j)
320: // * x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
321: // * end
322: // *
323: // * Define bounds on the components of x after j iterations of the loop:
324: // * M(j) = bound on x[1:j]
325: // * G(j) = bound on x[j+1:n]
326: // * Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
327: // *
328: // * Then for iteration j+1 we have
329: // * M(j+1) <= G(j) / | A(j+1,j+1) |
330: // * G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
331: // * <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
332: // *
333: // * where CNORM(j+1) is greater than or equal to the infinity-norm of
334: // * column j+1 of A, not counting the diagonal. Hence
335: // *
336: // * G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
337: // * 1<=i<=j
338: // * and
339: // *
340: // * |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
341: // * 1<=i< j
342: // *
343: // * Since |x(j)| <= M(j), we use the Level 2 BLAS routine DTRSV if the
344: // * reciprocal of the largest M(j), j=1,..,n, is larger than
345: // * max(underflow, 1/overflow).
346: // *
347: // * The bound on x(j) is also used to determine when a step in the
348: // * columnwise method can be performed without fear of overflow. If
349: // * the computed bound is greater than a large constant, x is scaled to
350: // * prevent overflow, but if the bound overflows, x is set to 0, x(j) to
351: // * 1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
352: // *
353: // * Similarly, a row-wise scheme is used to solve A'*x = b. The basic
354: // * algorithm for A upper triangular is
355: // *
356: // * for j = 1, ..., n
357: // * x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
358: // * end
359: // *
360: // * We simultaneously compute two bounds
361: // * G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
362: // * M(j) = bound on x(i), 1<=i<=j
363: // *
364: // * The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
365: // * add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
366: // * Then the bound on x(j) is
367: // *
368: // * M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
369: // *
370: // * <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
371: // * 1<=i<=j
372: // *
373: // * and we can safely call DTRSV if 1/M(n) and 1/G(n) are both greater
374: // * than max(underflow, 1/overflow).
375: // *
376: // * =====================================================================
377: // *
378: // * .. Parameters ..
379: // * ..
380: // * .. Local Scalars ..
381: // * ..
382: // * .. External Functions ..
383: // * ..
384: // * .. External Subroutines ..
385: // * ..
386: // * .. Intrinsic Functions ..
387: // INTRINSIC ABS, MAX, MIN;
388: // * ..
389: // * .. Executable Statements ..
390: // *
391:
392: #endregion
393:
394:
395: #region Body
396:
397: INFO = 0;
398: UPPER = this._lsame.Run(UPLO, "U");
399: NOTRAN = this._lsame.Run(TRANS, "N");
400: NOUNIT = this._lsame.Run(DIAG, "N");
401: // *
402: // * Test the input parameters.
403: // *
404: if (!UPPER && !this._lsame.Run(UPLO, "L"))
405: {
406: INFO = - 1;
407: }
408: else
409: {
410: if (!NOTRAN && !this._lsame.Run(TRANS, "T") && !this._lsame.Run(TRANS, "C"))
411: {
412: INFO = - 2;
413: }
414: else
415: {
416: if (!NOUNIT && !this._lsame.Run(DIAG, "U"))
417: {
418: INFO = - 3;
419: }
420: else
421: {
422: if (!this._lsame.Run(NORMIN, "Y") && !this._lsame.Run(NORMIN, "N"))
423: {
424: INFO = - 4;
425: }
426: else
427: {
428: if (N < 0)
429: {
430: INFO = - 5;
431: }
432: else
433: {
434: if (LDA < Math.Max(1, N))
435: {
436: INFO = - 7;
437: }
438: }
439: }
440: }
441: }
442: }
443: if (INFO != 0)
444: {
445: this._xerbla.Run("DLATRS", - INFO);
446: return;
447: }
448: // *
449: // * Quick return if possible
450: // *
451: if (N == 0) return;
452: // *
453: // * Determine machine dependent parameters to control overflow.
454: // *
455: SMLNUM = this._dlamch.Run("Safe minimum") / this._dlamch.Run("Precision");
456: BIGNUM = ONE / SMLNUM;
457: SCALE = ONE;
458: // *
459: if (this._lsame.Run(NORMIN, "N"))
460: {
461: // *
462: // * Compute the 1-norm of each column, not including the diagonal.
463: // *
464: if (UPPER)
465: {
466: // *
467: // * A is upper triangular.
468: // *
469: for (J = 1; J <= N; J++)
470: {
471: CNORM[J + o_cnorm] = this._dasum.Run(J - 1, A, 1+J * LDA + o_a, 1);
472: }
473: }
474: else
475: {
476: // *
477: // * A is lower triangular.
478: // *
479: for (J = 1; J <= N - 1; J++)
480: {
481: CNORM[J + o_cnorm] = this._dasum.Run(N - J, A, J + 1+J * LDA + o_a, 1);
482: }
483: CNORM[N + o_cnorm] = ZERO;
484: }
485: }
486: // *
487: // * Scale the column norms by TSCAL if the maximum element in CNORM is
488: // * greater than BIGNUM.
489: // *
490: IMAX = this._idamax.Run(N, CNORM, offset_cnorm, 1);
491: TMAX = CNORM[IMAX + o_cnorm];
492: if (TMAX <= BIGNUM)
493: {
494: TSCAL = ONE;
495: }
496: else
497: {
498: TSCAL = ONE / (SMLNUM * TMAX);
499: this._dscal.Run(N, TSCAL, ref CNORM, offset_cnorm, 1);
500: }
501: // *
502: // * Compute a bound on the computed solution vector to see if the
503: // * Level 2 BLAS routine DTRSV can be used.
504: // *
505: J = this._idamax.Run(N, X, offset_x, 1);
506: XMAX = Math.Abs(X[J + o_x]);
507: XBND = XMAX;
508: if (NOTRAN)
509: {
510: // *
511: // * Compute the growth in A * x = b.
512: // *
513: if (UPPER)
514: {
515: JFIRST = N;
516: JLAST = 1;
517: JINC = - 1;
518: }
519: else
520: {
521: JFIRST = 1;
522: JLAST = N;
523: JINC = 1;
524: }
525: // *
526: if (TSCAL != ONE)
527: {
528: GROW = ZERO;
529: goto LABEL50;
530: }
531: // *
532: if (NOUNIT)
533: {
534: // *
535: // * A is non-unit triangular.
536: // *
537: // * Compute GROW = 1/G(j) and XBND = 1/M(j).
538: // * Initially, G(0) = max{x(i), i=1,...,n}.
539: // *
540: GROW = ONE / Math.Max(XBND, SMLNUM);
541: XBND = GROW;
542: for (J = JFIRST; (JINC >= 0) ? (J <= JLAST) : (J >= JLAST); J += JINC)
543: {
544: // *
545: // * Exit the loop if the growth factor is too small.
546: // *
547: if (GROW <= SMLNUM) goto LABEL50;
548: // *
549: // * M(j) = G(j-1) / abs(A(j,j))
550: // *
551: TJJ = Math.Abs(A[J+J * LDA + o_a]);
552: XBND = Math.Min(XBND, Math.Min(ONE, TJJ) * GROW);
553: if (TJJ + CNORM[J + o_cnorm] >= SMLNUM)
554: {
555: // *
556: // * G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
557: // *
558: GROW = GROW * (TJJ / (TJJ + CNORM[J + o_cnorm]));
559: }
560: else
561: {
562: // *
563: // * G(j) could overflow, set GROW to 0.
564: // *
565: GROW = ZERO;
566: }
567: }
568: GROW = XBND;
569: }
570: else
571: {
572: // *
573: // * A is unit triangular.
574: // *
575: // * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
576: // *
577: GROW = Math.Min(ONE, ONE / Math.Max(XBND, SMLNUM));
578: for (J = JFIRST; (JINC >= 0) ? (J <= JLAST) : (J >= JLAST); J += JINC)
579: {
580: // *
581: // * Exit the loop if the growth factor is too small.
582: // *
583: if (GROW <= SMLNUM) goto LABEL50;
584: // *
585: // * G(j) = G(j-1)*( 1 + CNORM(j) )
586: // *
587: GROW = GROW * (ONE / (ONE + CNORM[J + o_cnorm]));
588: }
589: }
590: LABEL50:;
591: // *
592: }
593: else
594: {
595: // *
596: // * Compute the growth in A' * x = b.
597: // *
598: if (UPPER)
599: {
600: JFIRST = 1;
601: JLAST = N;
602: JINC = 1;
603: }
604: else
605: {
606: JFIRST = N;
607: JLAST = 1;
608: JINC = - 1;
609: }
610: // *
611: if (TSCAL != ONE)
612: {
613: GROW = ZERO;
614: goto LABEL80;
615: }
616: // *
617: if (NOUNIT)
618: {
619: // *
620: // * A is non-unit triangular.
621: // *
622: // * Compute GROW = 1/G(j) and XBND = 1/M(j).
623: // * Initially, M(0) = max{x(i), i=1,...,n}.
624: // *
625: GROW = ONE / Math.Max(XBND, SMLNUM);
626: XBND = GROW;
627: for (J = JFIRST; (JINC >= 0) ? (J <= JLAST) : (J >= JLAST); J += JINC)
628: {
629: // *
630: // * Exit the loop if the growth factor is too small.
631: // *
632: if (GROW <= SMLNUM) goto LABEL80;
633: // *
634: // * G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
635: // *
636: XJ = ONE + CNORM[J + o_cnorm];
637: GROW = Math.Min(GROW, XBND / XJ);
638: // *
639: // * M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
640: // *
641: TJJ = Math.Abs(A[J+J * LDA + o_a]);
642: if (XJ > TJJ) XBND = XBND * (TJJ / XJ);
643: }
644: GROW = Math.Min(GROW, XBND);
645: }
646: else
647: {
648: // *
649: // * A is unit triangular.
650: // *
651: // * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
652: // *
653: GROW = Math.Min(ONE, ONE / Math.Max(XBND, SMLNUM));
654: for (J = JFIRST; (JINC >= 0) ? (J <= JLAST) : (J >= JLAST); J += JINC)
655: {
656: // *
657: // * Exit the loop if the growth factor is too small.
658: // *
659: if (GROW <= SMLNUM) goto LABEL80;
660: // *
661: // * G(j) = ( 1 + CNORM(j) )*G(j-1)
662: // *
663: XJ = ONE + CNORM[J + o_cnorm];
664: GROW = GROW / XJ;
665: }
666: }
667: LABEL80:;
668: }
669: // *
670: if ((GROW * TSCAL) > SMLNUM)
671: {
672: // *
673: // * Use the Level 2 BLAS solve if the reciprocal of the bound on
674: // * elements of X is not too small.
675: // *
676: this._dtrsv.Run(UPLO, TRANS, DIAG, N, A, offset_a, LDA
677: , ref X, offset_x, 1);
678: }
679: else
680: {
681: // *
682: // * Use a Level 1 BLAS solve, scaling intermediate results.
683: // *
684: if (XMAX > BIGNUM)
685: {
686: // *
687: // * Scale X so that its components are less than or equal to
688: // * BIGNUM in absolute value.
689: // *
690: SCALE = BIGNUM / XMAX;
691: this._dscal.Run(N, SCALE, ref X, offset_x, 1);
692: XMAX = BIGNUM;
693: }
694: // *
695: if (NOTRAN)
696: {
697: // *
698: // * Solve A * x = b
699: // *
700: for (J = JFIRST; (JINC >= 0) ? (J <= JLAST) : (J >= JLAST); J += JINC)
701: {
702: // *
703: // * Compute x(j) = b(j) / A(j,j), scaling x if necessary.
704: // *
705: XJ = Math.Abs(X[J + o_x]);
706: if (NOUNIT)
707: {
708: TJJS = A[J+J * LDA + o_a] * TSCAL;
709: }
710: else
711: {
712: TJJS = TSCAL;
713: if (TSCAL == ONE) goto LABEL100;
714: }
715: TJJ = Math.Abs(TJJS);
716: if (TJJ > SMLNUM)
717: {
718: // *
719: // * abs(A(j,j)) > SMLNUM:
720: // *
721: if (TJJ < ONE)
722: {
723: if (XJ > TJJ * BIGNUM)
724: {
725: // *
726: // * Scale x by 1/b(j).
727: // *
728: REC = ONE / XJ;
729: this._dscal.Run(N, REC, ref X, offset_x, 1);
730: SCALE = SCALE * REC;
731: XMAX = XMAX * REC;
732: }
733: }
734: X[J + o_x] = X[J + o_x] / TJJS;
735: XJ = Math.Abs(X[J + o_x]);
736: }
737: else
738: {
739: if (TJJ > ZERO)
740: {
741: // *
742: // * 0 < abs(A(j,j)) <= SMLNUM:
743: // *
744: if (XJ > TJJ * BIGNUM)
745: {
746: // *
747: // * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
748: // * to avoid overflow when dividing by A(j,j).
749: // *
750: REC = (TJJ * BIGNUM) / XJ;
751: if (CNORM[J + o_cnorm] > ONE)
752: {
753: // *
754: // * Scale by 1/CNORM(j) to avoid overflow when
755: // * multiplying x(j) times column j.
756: // *
757: REC = REC / CNORM[J + o_cnorm];
758: }
759: this._dscal.Run(N, REC, ref X, offset_x, 1);
760: SCALE = SCALE * REC;
761: XMAX = XMAX * REC;
762: }
763: X[J + o_x] = X[J + o_x] / TJJS;
764: XJ = Math.Abs(X[J + o_x]);
765: }
766: else
767: {
768: // *
769: // * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
770: // * scale = 0, and compute a solution to A*x = 0.
771: // *
772: for (I = 1; I <= N; I++)
773: {
774: X[I + o_x] = ZERO;
775: }
776: X[J + o_x] = ONE;
777: XJ = ONE;
778: SCALE = ZERO;
779: XMAX = ZERO;
780: }
781: }
782: LABEL100:;
783: // *
784: // * Scale x if necessary to avoid overflow when adding a
785: // * multiple of column j of A.
786: // *
787: if (XJ > ONE)
788: {
789: REC = ONE / XJ;
790: if (CNORM[J + o_cnorm] > (BIGNUM - XMAX) * REC)
791: {
792: // *
793: // * Scale x by 1/(2*abs(x(j))).
794: // *
795: REC = REC * HALF;
796: this._dscal.Run(N, REC, ref X, offset_x, 1);
797: SCALE = SCALE * REC;
798: }
799: }
800: else
801: {
802: if (XJ * CNORM[J + o_cnorm] > (BIGNUM - XMAX))
803: {
804: // *
805: // * Scale x by 1/2.
806: // *
807: this._dscal.Run(N, HALF, ref X, offset_x, 1);
808: SCALE = SCALE * HALF;
809: }
810: }
811: // *
812: if (UPPER)
813: {
814: if (J > 1)
815: {
816: // *
817: // * Compute the update
818: // * x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j)
819: // *
820: this._daxpy.Run(J - 1, - X[J + o_x] * TSCAL, A, 1+J * LDA + o_a, 1, ref X, offset_x, 1);
821: I = this._idamax.Run(J - 1, X, offset_x, 1);
822: XMAX = Math.Abs(X[I + o_x]);
823: }
824: }
825: else
826: {
827: if (J < N)
828: {
829: // *
830: // * Compute the update
831: // * x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j)
832: // *
833: this._daxpy.Run(N - J, - X[J + o_x] * TSCAL, A, J + 1+J * LDA + o_a, 1, ref X, J + 1 + o_x, 1);
834: I = J + this._idamax.Run(N - J, X, J + 1 + o_x, 1);
835: XMAX = Math.Abs(X[I + o_x]);
836: }
837: }
838: }
839: // *
840: }
841: else
842: {
843: // *
844: // * Solve A' * x = b
845: // *
846: for (J = JFIRST; (JINC >= 0) ? (J <= JLAST) : (J >= JLAST); J += JINC)
847: {
848: // *
849: // * Compute x(j) = b(j) - sum A(k,j)*x(k).
850: // * k<>j
851: // *
852: XJ = Math.Abs(X[J + o_x]);
853: USCAL = TSCAL;
854: REC = ONE / Math.Max(XMAX, ONE);
855: if (CNORM[J + o_cnorm] > (BIGNUM - XJ) * REC)
856: {
857: // *
858: // * If x(j) could overflow, scale x by 1/(2*XMAX).
859: // *
860: REC = REC * HALF;
861: if (NOUNIT)
862: {
863: TJJS = A[J+J * LDA + o_a] * TSCAL;
864: }
865: else
866: {
867: TJJS = TSCAL;
868: }
869: TJJ = Math.Abs(TJJS);
870: if (TJJ > ONE)
871: {
872: // *
873: // * Divide by A(j,j) when scaling x if A(j,j) > 1.
874: // *
875: REC = Math.Min(ONE, REC * TJJ);
876: USCAL = USCAL / TJJS;
877: }
878: if (REC < ONE)
879: {
880: this._dscal.Run(N, REC, ref X, offset_x, 1);
881: SCALE = SCALE * REC;
882: XMAX = XMAX * REC;
883: }
884: }
885: // *
886: SUMJ = ZERO;
887: if (USCAL == ONE)
888: {
889: // *
890: // * If the scaling needed for A in the dot product is 1,
891: // * call DDOT to perform the dot product.
892: // *
893: if (UPPER)
894: {
895: SUMJ = this._ddot.Run(J - 1, A, 1+J * LDA + o_a, 1, X, offset_x, 1);
896: }
897: else
898: {
899: if (J < N)
900: {
901: SUMJ = this._ddot.Run(N - J, A, J + 1+J * LDA + o_a, 1, X, J + 1 + o_x, 1);
902: }
903: }
904: }
905: else
906: {
907: // *
908: // * Otherwise, use in-line code for the dot product.
909: // *
910: if (UPPER)
911: {
912: for (I = 1; I <= J - 1; I++)
913: {
914: SUMJ = SUMJ + (A[I+J * LDA + o_a] * USCAL) * X[I + o_x];
915: }
916: }
917: else
918: {
919: if (J < N)
920: {
921: for (I = J + 1; I <= N; I++)
922: {
923: SUMJ = SUMJ + (A[I+J * LDA + o_a] * USCAL) * X[I + o_x];
924: }
925: }
926: }
927: }
928: // *
929: if (USCAL == TSCAL)
930: {
931: // *
932: // * Compute x(j) := ( x(j) - sumj ) / A(j,j) if 1/A(j,j)
933: // * was not used to scale the dotproduct.
934: // *
935: X[J + o_x] = X[J + o_x] - SUMJ;
936: XJ = Math.Abs(X[J + o_x]);
937: if (NOUNIT)
938: {
939: TJJS = A[J+J * LDA + o_a] * TSCAL;
940: }
941: else
942: {
943: TJJS = TSCAL;
944: if (TSCAL == ONE) goto LABEL150;
945: }
946: // *
947: // * Compute x(j) = x(j) / A(j,j), scaling if necessary.
948: // *
949: TJJ = Math.Abs(TJJS);
950: if (TJJ > SMLNUM)
951: {
952: // *
953: // * abs(A(j,j)) > SMLNUM:
954: // *
955: if (TJJ < ONE)
956: {
957: if (XJ > TJJ * BIGNUM)
958: {
959: // *
960: // * Scale X by 1/abs(x(j)).
961: // *
962: REC = ONE / XJ;
963: this._dscal.Run(N, REC, ref X, offset_x, 1);
964: SCALE = SCALE * REC;
965: XMAX = XMAX * REC;
966: }
967: }
968: X[J + o_x] = X[J + o_x] / TJJS;
969: }
970: else
971: {
972: if (TJJ > ZERO)
973: {
974: // *
975: // * 0 < abs(A(j,j)) <= SMLNUM:
976: // *
977: if (XJ > TJJ * BIGNUM)
978: {
979: // *
980: // * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
981: // *
982: REC = (TJJ * BIGNUM) / XJ;
983: this._dscal.Run(N, REC, ref X, offset_x, 1);
984: SCALE = SCALE * REC;
985: XMAX = XMAX * REC;
986: }
987: X[J + o_x] = X[J + o_x] / TJJS;
988: }
989: else
990: {
991: // *
992: // * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
993: // * scale = 0, and compute a solution to A'*x = 0.
994: // *
995: for (I = 1; I <= N; I++)
996: {
997: X[I + o_x] = ZERO;
998: }
999: X[J + o_x] = ONE;
1000: SCALE = ZERO;
1001: XMAX = ZERO;
1002: }
1003: }
1004: LABEL150:;
1005: }
1006: else
1007: {
1008: // *
1009: // * Compute x(j) := x(j) / A(j,j) - sumj if the dot
1010: // * product has already been divided by 1/A(j,j).
1011: // *
1012: X[J + o_x] = X[J + o_x] / TJJS - SUMJ;
1013: }
1014: XMAX = Math.Max(XMAX, Math.Abs(X[J + o_x]));
1015: }
1016: }
1017: SCALE = SCALE / TSCAL;
1018: }
1019: // *
1020: // * Scale the column norms by 1/TSCAL for return.
1021: // *
1022: if (TSCAL != ONE)
1023: {
1024: this._dscal.Run(N, ONE / TSCAL, ref CNORM, offset_cnorm, 1);
1025: }
1026: // *
1027: return;
1028: // *
1029: // * End of DLATRS
1030: // *
1031:
1032: #endregion
1033:
1034: }
1035: }
1036: }