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: /// Purpose
21: /// =======
22: ///
23: /// DTRSM solves one of the matrix equations
24: ///
25: /// op( A )*X = alpha*B, or X*op( A ) = alpha*B,
26: ///
27: /// where alpha is a scalar, X and B are m by n matrices, A is a unit, or
28: /// non-unit, upper or lower triangular matrix and op( A ) is one of
29: ///
30: /// op( A ) = A or op( A ) = A'.
31: ///
32: /// The matrix X is overwritten on B.
33: ///
34: ///</summary>
35: public class DTRSM
36: {
37:
38:
39: #region Dependencies
40:
41: LSAME _lsame; XERBLA _xerbla;
42:
43: #endregion
44:
45:
46: #region Fields
47:
48: double TEMP = 0; int I = 0; int INFO = 0; int J = 0; int K = 0; int NROWA = 0; bool LSIDE = false; bool NOUNIT = false;
49: bool UPPER = false;const double ONE = 1.0E+0; const double ZERO = 0.0E+0;
50:
51: #endregion
52:
53: public DTRSM(LSAME lsame, XERBLA xerbla)
54: {
55:
56:
57: #region Set Dependencies
58:
59: this._lsame = lsame; this._xerbla = xerbla;
60:
61: #endregion
62:
63: }
64:
65: public DTRSM()
66: {
67:
68:
69: #region Dependencies (Initialization)
70:
71: LSAME lsame = new LSAME();
72: XERBLA xerbla = new XERBLA();
73:
74: #endregion
75:
76:
77: #region Set Dependencies
78:
79: this._lsame = lsame; this._xerbla = xerbla;
80:
81: #endregion
82:
83: }
84: /// <summary>
85: /// Purpose
86: /// =======
87: ///
88: /// DTRSM solves one of the matrix equations
89: ///
90: /// op( A )*X = alpha*B, or X*op( A ) = alpha*B,
91: ///
92: /// where alpha is a scalar, X and B are m by n matrices, A is a unit, or
93: /// non-unit, upper or lower triangular matrix and op( A ) is one of
94: ///
95: /// op( A ) = A or op( A ) = A'.
96: ///
97: /// The matrix X is overwritten on B.
98: ///
99: ///</summary>
100: /// <param name="SIDE">
101: /// - CHARACTER*1.
102: /// On entry, SIDE specifies whether op( A ) appears on the left
103: /// or right of X as follows:
104: ///
105: /// SIDE = 'L' or 'l' op( A )*X = alpha*B.
106: ///
107: /// SIDE = 'R' or 'r' X*op( A ) = alpha*B.
108: ///
109: /// Unchanged on exit.
110: ///</param>
111: /// <param name="UPLO">
112: /// - CHARACTER*1.
113: /// On entry, UPLO specifies whether the matrix A is an upper or
114: /// lower triangular matrix as follows:
115: ///
116: /// UPLO = 'U' or 'u' A is an upper triangular matrix.
117: ///
118: /// UPLO = 'L' or 'l' A is a lower triangular matrix.
119: ///
120: /// Unchanged on exit.
121: ///</param>
122: /// <param name="TRANSA">
123: /// - CHARACTER*1.
124: /// On entry, TRANSA specifies the form of op( A ) to be used in
125: /// the matrix multiplication as follows:
126: ///
127: /// TRANSA = 'N' or 'n' op( A ) = A.
128: ///
129: /// TRANSA = 'T' or 't' op( A ) = A'.
130: ///
131: /// TRANSA = 'C' or 'c' op( A ) = A'.
132: ///
133: /// Unchanged on exit.
134: ///</param>
135: /// <param name="DIAG">
136: /// - CHARACTER*1.
137: /// On entry, DIAG specifies whether or not A is unit triangular
138: /// as follows:
139: ///
140: /// DIAG = 'U' or 'u' A is assumed to be unit triangular.
141: ///
142: /// DIAG = 'N' or 'n' A is not assumed to be unit
143: /// triangular.
144: ///
145: /// Unchanged on exit.
146: ///</param>
147: /// <param name="M">
148: /// - INTEGER.
149: /// On entry, M specifies the number of rows of B. M must be at
150: /// least zero.
151: /// Unchanged on exit.
152: ///</param>
153: /// <param name="N">
154: /// - INTEGER.
155: /// On entry, N specifies the number of columns of B. N must be
156: /// at least zero.
157: /// Unchanged on exit.
158: ///</param>
159: /// <param name="ALPHA">
160: /// - DOUBLE PRECISION.
161: /// On entry, ALPHA specifies the scalar alpha. When alpha is
162: /// zero then A is not referenced and B need not be set before
163: /// entry.
164: /// Unchanged on exit.
165: ///</param>
166: /// <param name="A">
167: /// - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m
168: /// when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'.
169: /// Before entry with UPLO = 'U' or 'u', the leading k by k
170: /// upper triangular part of the array A must contain the upper
171: /// triangular matrix and the strictly lower triangular part of
172: /// A is not referenced.
173: /// Before entry with UPLO = 'L' or 'l', the leading k by k
174: /// lower triangular part of the array A must contain the lower
175: /// triangular matrix and the strictly upper triangular part of
176: /// A is not referenced.
177: /// Note that when DIAG = 'U' or 'u', the diagonal elements of
178: /// A are not referenced either, but are assumed to be unity.
179: /// Unchanged on exit.
180: ///</param>
181: /// <param name="LDA">
182: /// - INTEGER.
183: /// On entry, LDA specifies the first dimension of A as declared
184: /// in the calling (sub) program. When SIDE = 'L' or 'l' then
185: /// LDA must be at least max( 1, m ), when SIDE = 'R' or 'r'
186: /// then LDA must be at least max( 1, n ).
187: /// Unchanged on exit.
188: ///</param>
189: /// <param name="B">
190: /// - DOUBLE PRECISION array of DIMENSION ( LDB, n ).
191: /// Before entry, the leading m by n part of the array B must
192: /// contain the right-hand side matrix B, and on exit is
193: /// overwritten by the solution matrix X.
194: ///</param>
195: /// <param name="LDB">
196: /// - INTEGER.
197: /// On entry, LDB specifies the first dimension of B as declared
198: /// in the calling (sub) program. LDB must be at least
199: /// max( 1, m ).
200: /// Unchanged on exit.
201: ///
202: ///</param>
203: public void Run(string SIDE, string UPLO, string TRANSA, string DIAG, int M, int N
204: , double ALPHA, double[] A, int offset_a, int LDA, ref double[] B, int offset_b, int LDB)
205: {
206:
207: #region Array Index Correction
208:
209: int o_a = -1 - LDA + offset_a; int o_b = -1 - LDB + offset_b;
210:
211: #endregion
212:
213:
214: #region Strings
215:
216: SIDE = SIDE.Substring(0, 1); UPLO = UPLO.Substring(0, 1); TRANSA = TRANSA.Substring(0, 1);
217: DIAG = DIAG.Substring(0, 1);
218:
219: #endregion
220:
221:
222: #region Prolog
223:
224: // * .. Scalar Arguments ..
225: // * ..
226: // * .. Array Arguments ..
227: // * ..
228: // *
229: // * Purpose
230: // * =======
231: // *
232: // * DTRSM solves one of the matrix equations
233: // *
234: // * op( A )*X = alpha*B, or X*op( A ) = alpha*B,
235: // *
236: // * where alpha is a scalar, X and B are m by n matrices, A is a unit, or
237: // * non-unit, upper or lower triangular matrix and op( A ) is one of
238: // *
239: // * op( A ) = A or op( A ) = A'.
240: // *
241: // * The matrix X is overwritten on B.
242: // *
243: // * Arguments
244: // * ==========
245: // *
246: // * SIDE - CHARACTER*1.
247: // * On entry, SIDE specifies whether op( A ) appears on the left
248: // * or right of X as follows:
249: // *
250: // * SIDE = 'L' or 'l' op( A )*X = alpha*B.
251: // *
252: // * SIDE = 'R' or 'r' X*op( A ) = alpha*B.
253: // *
254: // * Unchanged on exit.
255: // *
256: // * UPLO - CHARACTER*1.
257: // * On entry, UPLO specifies whether the matrix A is an upper or
258: // * lower triangular matrix as follows:
259: // *
260: // * UPLO = 'U' or 'u' A is an upper triangular matrix.
261: // *
262: // * UPLO = 'L' or 'l' A is a lower triangular matrix.
263: // *
264: // * Unchanged on exit.
265: // *
266: // * TRANSA - CHARACTER*1.
267: // * On entry, TRANSA specifies the form of op( A ) to be used in
268: // * the matrix multiplication as follows:
269: // *
270: // * TRANSA = 'N' or 'n' op( A ) = A.
271: // *
272: // * TRANSA = 'T' or 't' op( A ) = A'.
273: // *
274: // * TRANSA = 'C' or 'c' op( A ) = A'.
275: // *
276: // * Unchanged on exit.
277: // *
278: // * DIAG - CHARACTER*1.
279: // * On entry, DIAG specifies whether or not A is unit triangular
280: // * as follows:
281: // *
282: // * DIAG = 'U' or 'u' A is assumed to be unit triangular.
283: // *
284: // * DIAG = 'N' or 'n' A is not assumed to be unit
285: // * triangular.
286: // *
287: // * Unchanged on exit.
288: // *
289: // * M - INTEGER.
290: // * On entry, M specifies the number of rows of B. M must be at
291: // * least zero.
292: // * Unchanged on exit.
293: // *
294: // * N - INTEGER.
295: // * On entry, N specifies the number of columns of B. N must be
296: // * at least zero.
297: // * Unchanged on exit.
298: // *
299: // * ALPHA - DOUBLE PRECISION.
300: // * On entry, ALPHA specifies the scalar alpha. When alpha is
301: // * zero then A is not referenced and B need not be set before
302: // * entry.
303: // * Unchanged on exit.
304: // *
305: // * A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m
306: // * when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'.
307: // * Before entry with UPLO = 'U' or 'u', the leading k by k
308: // * upper triangular part of the array A must contain the upper
309: // * triangular matrix and the strictly lower triangular part of
310: // * A is not referenced.
311: // * Before entry with UPLO = 'L' or 'l', the leading k by k
312: // * lower triangular part of the array A must contain the lower
313: // * triangular matrix and the strictly upper triangular part of
314: // * A is not referenced.
315: // * Note that when DIAG = 'U' or 'u', the diagonal elements of
316: // * A are not referenced either, but are assumed to be unity.
317: // * Unchanged on exit.
318: // *
319: // * LDA - INTEGER.
320: // * On entry, LDA specifies the first dimension of A as declared
321: // * in the calling (sub) program. When SIDE = 'L' or 'l' then
322: // * LDA must be at least max( 1, m ), when SIDE = 'R' or 'r'
323: // * then LDA must be at least max( 1, n ).
324: // * Unchanged on exit.
325: // *
326: // * B - DOUBLE PRECISION array of DIMENSION ( LDB, n ).
327: // * Before entry, the leading m by n part of the array B must
328: // * contain the right-hand side matrix B, and on exit is
329: // * overwritten by the solution matrix X.
330: // *
331: // * LDB - INTEGER.
332: // * On entry, LDB specifies the first dimension of B as declared
333: // * in the calling (sub) program. LDB must be at least
334: // * max( 1, m ).
335: // * Unchanged on exit.
336: // *
337: // *
338: // * Level 3 Blas routine.
339: // *
340: // *
341: // * -- Written on 8-February-1989.
342: // * Jack Dongarra, Argonne National Laboratory.
343: // * Iain Duff, AERE Harwell.
344: // * Jeremy Du Croz, Numerical Algorithms Group Ltd.
345: // * Sven Hammarling, Numerical Algorithms Group Ltd.
346: // *
347: // *
348: // * .. External Functions ..
349: // * ..
350: // * .. External Subroutines ..
351: // * ..
352: // * .. Intrinsic Functions ..
353: // INTRINSIC MAX;
354: // * ..
355: // * .. Local Scalars ..
356: // * ..
357: // * .. Parameters ..
358: // * ..
359: // *
360: // * Test the input parameters.
361: // *
362:
363: #endregion
364:
365:
366: #region Body
367:
368: LSIDE = this._lsame.Run(SIDE, "L");
369: if (LSIDE)
370: {
371: NROWA = M;
372: }
373: else
374: {
375: NROWA = N;
376: }
377: NOUNIT = this._lsame.Run(DIAG, "N");
378: UPPER = this._lsame.Run(UPLO, "U");
379: // *
380: INFO = 0;
381: if ((!LSIDE) && (!this._lsame.Run(SIDE, "R")))
382: {
383: INFO = 1;
384: }
385: else
386: {
387: if ((!UPPER) && (!this._lsame.Run(UPLO, "L")))
388: {
389: INFO = 2;
390: }
391: else
392: {
393: if ((!this._lsame.Run(TRANSA, "N")) && (!this._lsame.Run(TRANSA, "T")) && (!this._lsame.Run(TRANSA, "C")))
394: {
395: INFO = 3;
396: }
397: else
398: {
399: if ((!this._lsame.Run(DIAG, "U")) && (!this._lsame.Run(DIAG, "N")))
400: {
401: INFO = 4;
402: }
403: else
404: {
405: if (M < 0)
406: {
407: INFO = 5;
408: }
409: else
410: {
411: if (N < 0)
412: {
413: INFO = 6;
414: }
415: else
416: {
417: if (LDA < Math.Max(1, NROWA))
418: {
419: INFO = 9;
420: }
421: else
422: {
423: if (LDB < Math.Max(1, M))
424: {
425: INFO = 11;
426: }
427: }
428: }
429: }
430: }
431: }
432: }
433: }
434: if (INFO != 0)
435: {
436: this._xerbla.Run("DTRSM ", INFO);
437: return;
438: }
439: // *
440: // * Quick return if possible.
441: // *
442: if (N == 0) return;
443: // *
444: // * And when alpha.eq.zero.
445: // *
446: if (ALPHA == ZERO)
447: {
448: for (J = 1; J <= N; J++)
449: {
450: for (I = 1; I <= M; I++)
451: {
452: B[I+J * LDB + o_b] = ZERO;
453: }
454: }
455: return;
456: }
457: // *
458: // * Start the operations.
459: // *
460: if (LSIDE)
461: {
462: if (this._lsame.Run(TRANSA, "N"))
463: {
464: // *
465: // * Form B := alpha*inv( A )*B.
466: // *
467: if (UPPER)
468: {
469: for (J = 1; J <= N; J++)
470: {
471: if (ALPHA != ONE)
472: {
473: for (I = 1; I <= M; I++)
474: {
475: B[I+J * LDB + o_b] = ALPHA * B[I+J * LDB + o_b];
476: }
477: }
478: for (K = M; K >= 1; K += - 1)
479: {
480: if (B[K+J * LDB + o_b] != ZERO)
481: {
482: if (NOUNIT) B[K+J * LDB + o_b] = B[K+J * LDB + o_b] / A[K+K * LDA + o_a];
483: for (I = 1; I <= K - 1; I++)
484: {
485: B[I+J * LDB + o_b] = B[I+J * LDB + o_b] - B[K+J * LDB + o_b] * A[I+K * LDA + o_a];
486: }
487: }
488: }
489: }
490: }
491: else
492: {
493: for (J = 1; J <= N; J++)
494: {
495: if (ALPHA != ONE)
496: {
497: for (I = 1; I <= M; I++)
498: {
499: B[I+J * LDB + o_b] = ALPHA * B[I+J * LDB + o_b];
500: }
501: }
502: for (K = 1; K <= M; K++)
503: {
504: if (B[K+J * LDB + o_b] != ZERO)
505: {
506: if (NOUNIT) B[K+J * LDB + o_b] = B[K+J * LDB + o_b] / A[K+K * LDA + o_a];
507: for (I = K + 1; I <= M; I++)
508: {
509: B[I+J * LDB + o_b] = B[I+J * LDB + o_b] - B[K+J * LDB + o_b] * A[I+K * LDA + o_a];
510: }
511: }
512: }
513: }
514: }
515: }
516: else
517: {
518: // *
519: // * Form B := alpha*inv( A' )*B.
520: // *
521: if (UPPER)
522: {
523: for (J = 1; J <= N; J++)
524: {
525: for (I = 1; I <= M; I++)
526: {
527: TEMP = ALPHA * B[I+J * LDB + o_b];
528: for (K = 1; K <= I - 1; K++)
529: {
530: TEMP = TEMP - A[K+I * LDA + o_a] * B[K+J * LDB + o_b];
531: }
532: if (NOUNIT) TEMP = TEMP / A[I+I * LDA + o_a];
533: B[I+J * LDB + o_b] = TEMP;
534: }
535: }
536: }
537: else
538: {
539: for (J = 1; J <= N; J++)
540: {
541: for (I = M; I >= 1; I += - 1)
542: {
543: TEMP = ALPHA * B[I+J * LDB + o_b];
544: for (K = I + 1; K <= M; K++)
545: {
546: TEMP = TEMP - A[K+I * LDA + o_a] * B[K+J * LDB + o_b];
547: }
548: if (NOUNIT) TEMP = TEMP / A[I+I * LDA + o_a];
549: B[I+J * LDB + o_b] = TEMP;
550: }
551: }
552: }
553: }
554: }
555: else
556: {
557: if (this._lsame.Run(TRANSA, "N"))
558: {
559: // *
560: // * Form B := alpha*B*inv( A ).
561: // *
562: if (UPPER)
563: {
564: for (J = 1; J <= N; J++)
565: {
566: if (ALPHA != ONE)
567: {
568: for (I = 1; I <= M; I++)
569: {
570: B[I+J * LDB + o_b] = ALPHA * B[I+J * LDB + o_b];
571: }
572: }
573: for (K = 1; K <= J - 1; K++)
574: {
575: if (A[K+J * LDA + o_a] != ZERO)
576: {
577: for (I = 1; I <= M; I++)
578: {
579: B[I+J * LDB + o_b] = B[I+J * LDB + o_b] - A[K+J * LDA + o_a] * B[I+K * LDB + o_b];
580: }
581: }
582: }
583: if (NOUNIT)
584: {
585: TEMP = ONE / A[J+J * LDA + o_a];
586: for (I = 1; I <= M; I++)
587: {
588: B[I+J * LDB + o_b] = TEMP * B[I+J * LDB + o_b];
589: }
590: }
591: }
592: }
593: else
594: {
595: for (J = N; J >= 1; J += - 1)
596: {
597: if (ALPHA != ONE)
598: {
599: for (I = 1; I <= M; I++)
600: {
601: B[I+J * LDB + o_b] = ALPHA * B[I+J * LDB + o_b];
602: }
603: }
604: for (K = J + 1; K <= N; K++)
605: {
606: if (A[K+J * LDA + o_a] != ZERO)
607: {
608: for (I = 1; I <= M; I++)
609: {
610: B[I+J * LDB + o_b] = B[I+J * LDB + o_b] - A[K+J * LDA + o_a] * B[I+K * LDB + o_b];
611: }
612: }
613: }
614: if (NOUNIT)
615: {
616: TEMP = ONE / A[J+J * LDA + o_a];
617: for (I = 1; I <= M; I++)
618: {
619: B[I+J * LDB + o_b] = TEMP * B[I+J * LDB + o_b];
620: }
621: }
622: }
623: }
624: }
625: else
626: {
627: // *
628: // * Form B := alpha*B*inv( A' ).
629: // *
630: if (UPPER)
631: {
632: for (K = N; K >= 1; K += - 1)
633: {
634: if (NOUNIT)
635: {
636: TEMP = ONE / A[K+K * LDA + o_a];
637: for (I = 1; I <= M; I++)
638: {
639: B[I+K * LDB + o_b] = TEMP * B[I+K * LDB + o_b];
640: }
641: }
642: for (J = 1; J <= K - 1; J++)
643: {
644: if (A[J+K * LDA + o_a] != ZERO)
645: {
646: TEMP = A[J+K * LDA + o_a];
647: for (I = 1; I <= M; I++)
648: {
649: B[I+J * LDB + o_b] = B[I+J * LDB + o_b] - TEMP * B[I+K * LDB + o_b];
650: }
651: }
652: }
653: if (ALPHA != ONE)
654: {
655: for (I = 1; I <= M; I++)
656: {
657: B[I+K * LDB + o_b] = ALPHA * B[I+K * LDB + o_b];
658: }
659: }
660: }
661: }
662: else
663: {
664: for (K = 1; K <= N; K++)
665: {
666: if (NOUNIT)
667: {
668: TEMP = ONE / A[K+K * LDA + o_a];
669: for (I = 1; I <= M; I++)
670: {
671: B[I+K * LDB + o_b] = TEMP * B[I+K * LDB + o_b];
672: }
673: }
674: for (J = K + 1; J <= N; J++)
675: {
676: if (A[J+K * LDA + o_a] != ZERO)
677: {
678: TEMP = A[J+K * LDA + o_a];
679: for (I = 1; I <= M; I++)
680: {
681: B[I+J * LDB + o_b] = B[I+J * LDB + o_b] - TEMP * B[I+K * LDB + o_b];
682: }
683: }
684: }
685: if (ALPHA != ONE)
686: {
687: for (I = 1; I <= M; I++)
688: {
689: B[I+K * LDB + o_b] = ALPHA * B[I+K * LDB + o_b];
690: }
691: }
692: }
693: }
694: }
695: }
696: // *
697: return;
698: // *
699: // * End of DTRSM .
700: // *
701:
702: #endregion
703:
704: }
705: }
706: }