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: /// October 31, 1992
24: /// Purpose
25: /// =======
26: ///
27: /// DLANTR returns the value of the one norm, or the Frobenius norm, or
28: /// the infinity norm, or the element of largest absolute value of a
29: /// trapezoidal or triangular matrix A.
30: ///
31: /// Description
32: /// ===========
33: ///
34: /// DLANTR returns the value
35: ///
36: /// DLANTR = ( max(abs(A(i,j))), NORM = 'M' or 'm'
37: /// (
38: /// ( norm1(A), NORM = '1', 'O' or 'o'
39: /// (
40: /// ( normI(A), NORM = 'I' or 'i'
41: /// (
42: /// ( normF(A), NORM = 'F', 'f', 'E' or 'e'
43: ///
44: /// where norm1 denotes the one norm of a matrix (maximum column sum),
45: /// normI denotes the infinity norm of a matrix (maximum row sum) and
46: /// normF denotes the Frobenius norm of a matrix (square root of sum of
47: /// squares). Note that max(abs(A(i,j))) is not a matrix norm.
48: ///
49: ///</summary>
50: public class DLANTR
51: {
52:
53:
54: #region Dependencies
55:
56: DLASSQ _dlassq; LSAME _lsame;
57:
58: #endregion
59:
60:
61: #region Fields
62:
63: const double ONE = 1.0E+0; const double ZERO = 0.0E+0; bool UDIAG = false; int I = 0; int J = 0; double SCALE = 0;
64: double SUM = 0;double VALUE = 0;
65:
66: #endregion
67:
68: public DLANTR(DLASSQ dlassq, LSAME lsame)
69: {
70:
71:
72: #region Set Dependencies
73:
74: this._dlassq = dlassq; this._lsame = lsame;
75:
76: #endregion
77:
78: }
79:
80: public DLANTR()
81: {
82:
83:
84: #region Dependencies (Initialization)
85:
86: DLASSQ dlassq = new DLASSQ();
87: LSAME lsame = new LSAME();
88:
89: #endregion
90:
91:
92: #region Set Dependencies
93:
94: this._dlassq = dlassq; this._lsame = lsame;
95:
96: #endregion
97:
98: }
99: /// <summary>
100: /// Purpose
101: /// =======
102: ///
103: /// DLANTR returns the value of the one norm, or the Frobenius norm, or
104: /// the infinity norm, or the element of largest absolute value of a
105: /// trapezoidal or triangular matrix A.
106: ///
107: /// Description
108: /// ===========
109: ///
110: /// DLANTR returns the value
111: ///
112: /// DLANTR = ( max(abs(A(i,j))), NORM = 'M' or 'm'
113: /// (
114: /// ( norm1(A), NORM = '1', 'O' or 'o'
115: /// (
116: /// ( normI(A), NORM = 'I' or 'i'
117: /// (
118: /// ( normF(A), NORM = 'F', 'f', 'E' or 'e'
119: ///
120: /// where norm1 denotes the one norm of a matrix (maximum column sum),
121: /// normI denotes the infinity norm of a matrix (maximum row sum) and
122: /// normF denotes the Frobenius norm of a matrix (square root of sum of
123: /// squares). Note that max(abs(A(i,j))) is not a matrix norm.
124: ///
125: ///</summary>
126: /// <param name="NORM">
127: /// (input) CHARACTER*1
128: /// Specifies the value to be returned in DLANTR as described
129: /// above.
130: ///</param>
131: /// <param name="UPLO">
132: /// (input) CHARACTER*1
133: /// Specifies whether the matrix A is upper or lower trapezoidal.
134: /// = 'U': Upper trapezoidal
135: /// = 'L': Lower trapezoidal
136: /// Note that A is triangular instead of trapezoidal if M = N.
137: ///</param>
138: /// <param name="DIAG">
139: /// (input) CHARACTER*1
140: /// Specifies whether or not the matrix A has unit diagonal.
141: /// = 'N': Non-unit diagonal
142: /// = 'U': Unit diagonal
143: ///</param>
144: /// <param name="M">
145: /// (input) INTEGER
146: /// The number of rows of the matrix A. M .GE. 0, and if
147: /// UPLO = 'U', M .LE. N. When M = 0, DLANTR is set to zero.
148: ///</param>
149: /// <param name="N">
150: /// (input) INTEGER
151: /// The number of columns of the matrix A. N .GE. 0, and if
152: /// UPLO = 'L', N .LE. M. When N = 0, DLANTR is set to zero.
153: ///</param>
154: /// <param name="A">
155: /// (input) DOUBLE PRECISION array, dimension (LDA,N)
156: /// The trapezoidal matrix A (A is triangular if M = N).
157: /// If UPLO = 'U', the leading m by n upper trapezoidal part of
158: /// the array A contains the upper trapezoidal matrix, and the
159: /// strictly lower triangular part of A is not referenced.
160: /// If UPLO = 'L', the leading m by n lower trapezoidal part of
161: /// the array A contains the lower trapezoidal matrix, and the
162: /// strictly upper triangular part of A is not referenced. Note
163: /// that when DIAG = 'U', the diagonal elements of A are not
164: /// referenced and are assumed to be one.
165: ///</param>
166: /// <param name="LDA">
167: /// (input) INTEGER
168: /// The leading dimension of the array A. LDA .GE. max(M,1).
169: ///</param>
170: /// <param name="WORK">
171: /// (workspace) DOUBLE PRECISION array, dimension (LWORK),
172: /// where LWORK .GE. M when NORM = 'I'; otherwise, WORK is not
173: /// referenced.
174: ///</param>
175: public double Run(string NORM, string UPLO, string DIAG, int M, int N, double[] A, int offset_a
176: , int LDA, ref double[] WORK, int offset_work)
177: {
178: double dlantr = 0;
179:
180: #region Array Index Correction
181:
182: int o_a = -1 - LDA + offset_a; int o_work = -1 + offset_work;
183:
184: #endregion
185:
186:
187: #region Strings
188:
189: NORM = NORM.Substring(0, 1); UPLO = UPLO.Substring(0, 1); DIAG = DIAG.Substring(0, 1);
190:
191: #endregion
192:
193:
194: #region Prolog
195:
196: // *
197: // * -- LAPACK auxiliary routine (version 3.0) --
198: // * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
199: // * Courant Institute, Argonne National Lab, and Rice University
200: // * October 31, 1992
201: // *
202: // * .. Scalar Arguments ..
203: // * ..
204: // * .. Array Arguments ..
205: // * ..
206: // *
207: // * Purpose
208: // * =======
209: // *
210: // * DLANTR returns the value of the one norm, or the Frobenius norm, or
211: // * the infinity norm, or the element of largest absolute value of a
212: // * trapezoidal or triangular matrix A.
213: // *
214: // * Description
215: // * ===========
216: // *
217: // * DLANTR returns the value
218: // *
219: // * DLANTR = ( max(abs(A(i,j))), NORM = 'M' or 'm'
220: // * (
221: // * ( norm1(A), NORM = '1', 'O' or 'o'
222: // * (
223: // * ( normI(A), NORM = 'I' or 'i'
224: // * (
225: // * ( normF(A), NORM = 'F', 'f', 'E' or 'e'
226: // *
227: // * where norm1 denotes the one norm of a matrix (maximum column sum),
228: // * normI denotes the infinity norm of a matrix (maximum row sum) and
229: // * normF denotes the Frobenius norm of a matrix (square root of sum of
230: // * squares). Note that max(abs(A(i,j))) is not a matrix norm.
231: // *
232: // * Arguments
233: // * =========
234: // *
235: // * NORM (input) CHARACTER*1
236: // * Specifies the value to be returned in DLANTR as described
237: // * above.
238: // *
239: // * UPLO (input) CHARACTER*1
240: // * Specifies whether the matrix A is upper or lower trapezoidal.
241: // * = 'U': Upper trapezoidal
242: // * = 'L': Lower trapezoidal
243: // * Note that A is triangular instead of trapezoidal if M = N.
244: // *
245: // * DIAG (input) CHARACTER*1
246: // * Specifies whether or not the matrix A has unit diagonal.
247: // * = 'N': Non-unit diagonal
248: // * = 'U': Unit diagonal
249: // *
250: // * M (input) INTEGER
251: // * The number of rows of the matrix A. M >= 0, and if
252: // * UPLO = 'U', M <= N. When M = 0, DLANTR is set to zero.
253: // *
254: // * N (input) INTEGER
255: // * The number of columns of the matrix A. N >= 0, and if
256: // * UPLO = 'L', N <= M. When N = 0, DLANTR is set to zero.
257: // *
258: // * A (input) DOUBLE PRECISION array, dimension (LDA,N)
259: // * The trapezoidal matrix A (A is triangular if M = N).
260: // * If UPLO = 'U', the leading m by n upper trapezoidal part of
261: // * the array A contains the upper trapezoidal matrix, and the
262: // * strictly lower triangular part of A is not referenced.
263: // * If UPLO = 'L', the leading m by n lower trapezoidal part of
264: // * the array A contains the lower trapezoidal matrix, and the
265: // * strictly upper triangular part of A is not referenced. Note
266: // * that when DIAG = 'U', the diagonal elements of A are not
267: // * referenced and are assumed to be one.
268: // *
269: // * LDA (input) INTEGER
270: // * The leading dimension of the array A. LDA >= max(M,1).
271: // *
272: // * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK),
273: // * where LWORK >= M when NORM = 'I'; otherwise, WORK is not
274: // * referenced.
275: // *
276: // * =====================================================================
277: // *
278: // * .. Parameters ..
279: // * ..
280: // * .. Local Scalars ..
281: // * ..
282: // * .. External Subroutines ..
283: // * ..
284: // * .. External Functions ..
285: // * ..
286: // * .. Intrinsic Functions ..
287: // INTRINSIC ABS, MAX, MIN, SQRT;
288: // * ..
289: // * .. Executable Statements ..
290: // *
291:
292: #endregion
293:
294:
295: #region Body
296:
297: if (Math.Min(M, N) == 0)
298: {
299: VALUE = ZERO;
300: }
301: else
302: {
303: if (this._lsame.Run(NORM, "M"))
304: {
305: // *
306: // * Find max(abs(A(i,j))).
307: // *
308: if (this._lsame.Run(DIAG, "U"))
309: {
310: VALUE = ONE;
311: if (this._lsame.Run(UPLO, "U"))
312: {
313: for (J = 1; J <= N; J++)
314: {
315: for (I = 1; I <= Math.Min(M, J - 1); I++)
316: {
317: VALUE = Math.Max(VALUE, Math.Abs(A[I+J * LDA + o_a]));
318: }
319: }
320: }
321: else
322: {
323: for (J = 1; J <= N; J++)
324: {
325: for (I = J + 1; I <= M; I++)
326: {
327: VALUE = Math.Max(VALUE, Math.Abs(A[I+J * LDA + o_a]));
328: }
329: }
330: }
331: }
332: else
333: {
334: VALUE = ZERO;
335: if (this._lsame.Run(UPLO, "U"))
336: {
337: for (J = 1; J <= N; J++)
338: {
339: for (I = 1; I <= Math.Min(M, J); I++)
340: {
341: VALUE = Math.Max(VALUE, Math.Abs(A[I+J * LDA + o_a]));
342: }
343: }
344: }
345: else
346: {
347: for (J = 1; J <= N; J++)
348: {
349: for (I = J; I <= M; I++)
350: {
351: VALUE = Math.Max(VALUE, Math.Abs(A[I+J * LDA + o_a]));
352: }
353: }
354: }
355: }
356: }
357: else
358: {
359: if ((this._lsame.Run(NORM, "O")) || (NORM == "1"))
360: {
361: // *
362: // * Find norm1(A).
363: // *
364: VALUE = ZERO;
365: UDIAG = this._lsame.Run(DIAG, "U");
366: if (this._lsame.Run(UPLO, "U"))
367: {
368: for (J = 1; J <= N; J++)
369: {
370: if ((UDIAG) && (J <= M))
371: {
372: SUM = ONE;
373: for (I = 1; I <= J - 1; I++)
374: {
375: SUM = SUM + Math.Abs(A[I+J * LDA + o_a]);
376: }
377: }
378: else
379: {
380: SUM = ZERO;
381: for (I = 1; I <= Math.Min(M, J); I++)
382: {
383: SUM = SUM + Math.Abs(A[I+J * LDA + o_a]);
384: }
385: }
386: VALUE = Math.Max(VALUE, SUM);
387: }
388: }
389: else
390: {
391: for (J = 1; J <= N; J++)
392: {
393: if (UDIAG)
394: {
395: SUM = ONE;
396: for (I = J + 1; I <= M; I++)
397: {
398: SUM = SUM + Math.Abs(A[I+J * LDA + o_a]);
399: }
400: }
401: else
402: {
403: SUM = ZERO;
404: for (I = J; I <= M; I++)
405: {
406: SUM = SUM + Math.Abs(A[I+J * LDA + o_a]);
407: }
408: }
409: VALUE = Math.Max(VALUE, SUM);
410: }
411: }
412: }
413: else
414: {
415: if (this._lsame.Run(NORM, "I"))
416: {
417: // *
418: // * Find normI(A).
419: // *
420: if (this._lsame.Run(UPLO, "U"))
421: {
422: if (this._lsame.Run(DIAG, "U"))
423: {
424: for (I = 1; I <= M; I++)
425: {
426: WORK[I + o_work] = ONE;
427: }
428: for (J = 1; J <= N; J++)
429: {
430: for (I = 1; I <= Math.Min(M, J - 1); I++)
431: {
432: WORK[I + o_work] = WORK[I + o_work] + Math.Abs(A[I+J * LDA + o_a]);
433: }
434: }
435: }
436: else
437: {
438: for (I = 1; I <= M; I++)
439: {
440: WORK[I + o_work] = ZERO;
441: }
442: for (J = 1; J <= N; J++)
443: {
444: for (I = 1; I <= Math.Min(M, J); I++)
445: {
446: WORK[I + o_work] = WORK[I + o_work] + Math.Abs(A[I+J * LDA + o_a]);
447: }
448: }
449: }
450: }
451: else
452: {
453: if (this._lsame.Run(DIAG, "U"))
454: {
455: for (I = 1; I <= N; I++)
456: {
457: WORK[I + o_work] = ONE;
458: }
459: for (I = N + 1; I <= M; I++)
460: {
461: WORK[I + o_work] = ZERO;
462: }
463: for (J = 1; J <= N; J++)
464: {
465: for (I = J + 1; I <= M; I++)
466: {
467: WORK[I + o_work] = WORK[I + o_work] + Math.Abs(A[I+J * LDA + o_a]);
468: }
469: }
470: }
471: else
472: {
473: for (I = 1; I <= M; I++)
474: {
475: WORK[I + o_work] = ZERO;
476: }
477: for (J = 1; J <= N; J++)
478: {
479: for (I = J; I <= M; I++)
480: {
481: WORK[I + o_work] = WORK[I + o_work] + Math.Abs(A[I+J * LDA + o_a]);
482: }
483: }
484: }
485: }
486: VALUE = ZERO;
487: for (I = 1; I <= M; I++)
488: {
489: VALUE = Math.Max(VALUE, WORK[I + o_work]);
490: }
491: }
492: else
493: {
494: if ((this._lsame.Run(NORM, "F")) || (this._lsame.Run(NORM, "E")))
495: {
496: // *
497: // * Find normF(A).
498: // *
499: if (this._lsame.Run(UPLO, "U"))
500: {
501: if (this._lsame.Run(DIAG, "U"))
502: {
503: SCALE = ONE;
504: SUM = Math.Min(M, N);
505: for (J = 2; J <= N; J++)
506: {
507: this._dlassq.Run(Math.Min(M, J - 1), A, 1+J * LDA + o_a, 1, ref SCALE, ref SUM);
508: }
509: }
510: else
511: {
512: SCALE = ZERO;
513: SUM = ONE;
514: for (J = 1; J <= N; J++)
515: {
516: this._dlassq.Run(Math.Min(M, J), A, 1+J * LDA + o_a, 1, ref SCALE, ref SUM);
517: }
518: }
519: }
520: else
521: {
522: if (this._lsame.Run(DIAG, "U"))
523: {
524: SCALE = ONE;
525: SUM = Math.Min(M, N);
526: for (J = 1; J <= N; J++)
527: {
528: this._dlassq.Run(M - J, A, Math.Min(M, J + 1)+J * LDA + o_a, 1, ref SCALE, ref SUM);
529: }
530: }
531: else
532: {
533: SCALE = ZERO;
534: SUM = ONE;
535: for (J = 1; J <= N; J++)
536: {
537: this._dlassq.Run(M - J + 1, A, J+J * LDA + o_a, 1, ref SCALE, ref SUM);
538: }
539: }
540: }
541: VALUE = SCALE * Math.Sqrt(SUM);
542: }
543: }
544: }
545: }
546: }
547: // *
548: dlantr = VALUE;
549: return dlantr;
550: // *
551: // * End of DLANTR
552: // *
553:
554: #endregion
555:
556: }
557: }
558: }