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: /// DLARFB applies a real block reflector H or its transpose H' to a
27: /// real m by n matrix C, from either the left or the right.
28: ///
29: ///</summary>
30: public class DLARFB
31: {
32:
33:
34: #region Dependencies
35:
36: LSAME _lsame; DCOPY _dcopy; DGEMM _dgemm; DTRMM _dtrmm;
37:
38: #endregion
39:
40:
41: #region Fields
42:
43: const double ONE = 1.0E+0; string TRANST = new string(' ', 1); int I = 0; int J = 0;
44:
45: #endregion
46:
47: public DLARFB(LSAME lsame, DCOPY dcopy, DGEMM dgemm, DTRMM dtrmm)
48: {
49:
50:
51: #region Set Dependencies
52:
53: this._lsame = lsame; this._dcopy = dcopy; this._dgemm = dgemm; this._dtrmm = dtrmm;
54:
55: #endregion
56:
57: }
58:
59: public DLARFB()
60: {
61:
62:
63: #region Dependencies (Initialization)
64:
65: LSAME lsame = new LSAME();
66: DCOPY dcopy = new DCOPY();
67: XERBLA xerbla = new XERBLA();
68: DGEMM dgemm = new DGEMM(lsame, xerbla);
69: DTRMM dtrmm = new DTRMM(lsame, xerbla);
70:
71: #endregion
72:
73:
74: #region Set Dependencies
75:
76: this._lsame = lsame; this._dcopy = dcopy; this._dgemm = dgemm; this._dtrmm = dtrmm;
77:
78: #endregion
79:
80: }
81: /// <summary>
82: /// Purpose
83: /// =======
84: ///
85: /// DLARFB applies a real block reflector H or its transpose H' to a
86: /// real m by n matrix C, from either the left or the right.
87: ///
88: ///</summary>
89: /// <param name="SIDE">
90: /// (input) CHARACTER*1
91: /// = 'L': apply H or H' from the Left
92: /// = 'R': apply H or H' from the Right
93: ///</param>
94: /// <param name="TRANS">
95: /// (input) CHARACTER*1
96: /// = 'N': apply H (No transpose)
97: /// = 'T': apply H' (Transpose)
98: ///</param>
99: /// <param name="DIRECT">
100: /// (input) CHARACTER*1
101: /// Indicates how H is formed from a product of elementary
102: /// reflectors
103: /// = 'F': H = H(1) H(2) . . . H(k) (Forward)
104: /// = 'B': H = H(k) . . . H(2) H(1) (Backward)
105: ///</param>
106: /// <param name="STOREV">
107: /// (input) CHARACTER*1
108: /// Indicates how the vectors which define the elementary
109: /// reflectors are stored:
110: /// = 'C': Columnwise
111: /// = 'R': Rowwise
112: ///</param>
113: /// <param name="M">
114: /// (input) INTEGER
115: /// The number of rows of the matrix C.
116: ///</param>
117: /// <param name="N">
118: /// (input) INTEGER
119: /// The number of columns of the matrix C.
120: ///</param>
121: /// <param name="K">
122: /// (input) INTEGER
123: /// The order of the matrix T (= the number of elementary
124: /// reflectors whose product defines the block reflector).
125: ///</param>
126: /// <param name="V">
127: /// (input) DOUBLE PRECISION array, dimension
128: /// (LDV,K) if STOREV = 'C'
129: /// (LDV,M) if STOREV = 'R' and SIDE = 'L'
130: /// (LDV,N) if STOREV = 'R' and SIDE = 'R'
131: /// The matrix V. See further details.
132: ///</param>
133: /// <param name="LDV">
134: /// (input) INTEGER
135: /// The leading dimension of the array V.
136: /// If STOREV = 'C' and SIDE = 'L', LDV .GE. max(1,M);
137: /// if STOREV = 'C' and SIDE = 'R', LDV .GE. max(1,N);
138: /// if STOREV = 'R', LDV .GE. K.
139: ///</param>
140: /// <param name="T">
141: /// (input) DOUBLE PRECISION array, dimension (LDT,K)
142: /// The triangular k by k matrix T in the representation of the
143: /// block reflector.
144: ///</param>
145: /// <param name="LDT">
146: /// (input) INTEGER
147: /// The leading dimension of the array T. LDT .GE. K.
148: ///</param>
149: /// <param name="C">
150: /// (input/output) DOUBLE PRECISION array, dimension (LDC,N)
151: /// On entry, the m by n matrix C.
152: /// On exit, C is overwritten by H*C or H'*C or C*H or C*H'.
153: ///</param>
154: /// <param name="LDC">
155: /// (input) INTEGER
156: /// The leading dimension of the array C. LDA .GE. max(1,M).
157: ///</param>
158: /// <param name="WORK">
159: /// (workspace) DOUBLE PRECISION array, dimension (LDWORK,K)
160: ///</param>
161: /// <param name="LDWORK">
162: /// (input) INTEGER
163: /// The leading dimension of the array WORK.
164: /// If SIDE = 'L', LDWORK .GE. max(1,N);
165: /// if SIDE = 'R', LDWORK .GE. max(1,M).
166: ///</param>
167: public void Run(string SIDE, string TRANS, string DIRECT, string STOREV, int M, int N
168: , int K, double[] V, int offset_v, int LDV, double[] T, int offset_t, int LDT, ref double[] C, int offset_c
169: , int LDC, ref double[] WORK, int offset_work, int LDWORK)
170: {
171:
172: #region Array Index Correction
173:
174: int o_v = -1 - LDV + offset_v; int o_t = -1 - LDT + offset_t; int o_c = -1 - LDC + offset_c;
175: int o_work = -1 - LDWORK + offset_work;
176:
177: #endregion
178:
179:
180: #region Strings
181:
182: SIDE = SIDE.Substring(0, 1); TRANS = TRANS.Substring(0, 1); DIRECT = DIRECT.Substring(0, 1);
183: STOREV = STOREV.Substring(0, 1);
184:
185: #endregion
186:
187:
188: #region Prolog
189:
190: // *
191: // * -- LAPACK auxiliary routine (version 3.1) --
192: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
193: // * November 2006
194: // *
195: // * .. Scalar Arguments ..
196: // * ..
197: // * .. Array Arguments ..
198: // * ..
199: // *
200: // * Purpose
201: // * =======
202: // *
203: // * DLARFB applies a real block reflector H or its transpose H' to a
204: // * real m by n matrix C, from either the left or the right.
205: // *
206: // * Arguments
207: // * =========
208: // *
209: // * SIDE (input) CHARACTER*1
210: // * = 'L': apply H or H' from the Left
211: // * = 'R': apply H or H' from the Right
212: // *
213: // * TRANS (input) CHARACTER*1
214: // * = 'N': apply H (No transpose)
215: // * = 'T': apply H' (Transpose)
216: // *
217: // * DIRECT (input) CHARACTER*1
218: // * Indicates how H is formed from a product of elementary
219: // * reflectors
220: // * = 'F': H = H(1) H(2) . . . H(k) (Forward)
221: // * = 'B': H = H(k) . . . H(2) H(1) (Backward)
222: // *
223: // * STOREV (input) CHARACTER*1
224: // * Indicates how the vectors which define the elementary
225: // * reflectors are stored:
226: // * = 'C': Columnwise
227: // * = 'R': Rowwise
228: // *
229: // * M (input) INTEGER
230: // * The number of rows of the matrix C.
231: // *
232: // * N (input) INTEGER
233: // * The number of columns of the matrix C.
234: // *
235: // * K (input) INTEGER
236: // * The order of the matrix T (= the number of elementary
237: // * reflectors whose product defines the block reflector).
238: // *
239: // * V (input) DOUBLE PRECISION array, dimension
240: // * (LDV,K) if STOREV = 'C'
241: // * (LDV,M) if STOREV = 'R' and SIDE = 'L'
242: // * (LDV,N) if STOREV = 'R' and SIDE = 'R'
243: // * The matrix V. See further details.
244: // *
245: // * LDV (input) INTEGER
246: // * The leading dimension of the array V.
247: // * If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
248: // * if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
249: // * if STOREV = 'R', LDV >= K.
250: // *
251: // * T (input) DOUBLE PRECISION array, dimension (LDT,K)
252: // * The triangular k by k matrix T in the representation of the
253: // * block reflector.
254: // *
255: // * LDT (input) INTEGER
256: // * The leading dimension of the array T. LDT >= K.
257: // *
258: // * C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
259: // * On entry, the m by n matrix C.
260: // * On exit, C is overwritten by H*C or H'*C or C*H or C*H'.
261: // *
262: // * LDC (input) INTEGER
263: // * The leading dimension of the array C. LDA >= max(1,M).
264: // *
265: // * WORK (workspace) DOUBLE PRECISION array, dimension (LDWORK,K)
266: // *
267: // * LDWORK (input) INTEGER
268: // * The leading dimension of the array WORK.
269: // * If SIDE = 'L', LDWORK >= max(1,N);
270: // * if SIDE = 'R', LDWORK >= max(1,M).
271: // *
272: // * =====================================================================
273: // *
274: // * .. Parameters ..
275: // * ..
276: // * .. Local Scalars ..
277: // * ..
278: // * .. External Functions ..
279: // * ..
280: // * .. External Subroutines ..
281: // * ..
282: // * .. Executable Statements ..
283: // *
284: // * Quick return if possible
285: // *
286:
287: #endregion
288:
289:
290: #region Body
291:
292: if (M <= 0 || N <= 0) return;
293: // *
294: if (this._lsame.Run(TRANS, "N"))
295: {
296: FortranLib.Copy(ref TRANST , "T");
297: }
298: else
299: {
300: FortranLib.Copy(ref TRANST , "N");
301: }
302: // *
303: if (this._lsame.Run(STOREV, "C"))
304: {
305: // *
306: if (this._lsame.Run(DIRECT, "F"))
307: {
308: // *
309: // * Let V = ( V1 ) (first K rows)
310: // * ( V2 )
311: // * where V1 is unit lower triangular.
312: // *
313: if (this._lsame.Run(SIDE, "L"))
314: {
315: // *
316: // * Form H * C or H' * C where C = ( C1 )
317: // * ( C2 )
318: // *
319: // * W := C' * V = (C1'*V1 + C2'*V2) (stored in WORK)
320: // *
321: // * W := C1'
322: // *
323: for (J = 1; J <= K; J++)
324: {
325: this._dcopy.Run(N, C, J+1 * LDC + o_c, LDC, ref WORK, 1+J * LDWORK + o_work, 1);
326: }
327: // *
328: // * W := W * V1
329: // *
330: this._dtrmm.Run("Right", "Lower", "No transpose", "Unit", N, K
331: , ONE, V, offset_v, LDV, ref WORK, offset_work, LDWORK);
332: if (M > K)
333: {
334: // *
335: // * W := W + C2'*V2
336: // *
337: this._dgemm.Run("Transpose", "No transpose", N, K, M - K, ONE
338: , C, K + 1+1 * LDC + o_c, LDC, V, K + 1+1 * LDV + o_v, LDV, ONE, ref WORK, offset_work
339: , LDWORK);
340: }
341: // *
342: // * W := W * T' or W * T
343: // *
344: this._dtrmm.Run("Right", "Upper", TRANST, "Non-unit", N, K
345: , ONE, T, offset_t, LDT, ref WORK, offset_work, LDWORK);
346: // *
347: // * C := C - V * W'
348: // *
349: if (M > K)
350: {
351: // *
352: // * C2 := C2 - V2 * W'
353: // *
354: this._dgemm.Run("No transpose", "Transpose", M - K, N, K, - ONE
355: , V, K + 1+1 * LDV + o_v, LDV, WORK, offset_work, LDWORK, ONE, ref C, K + 1+1 * LDC + o_c
356: , LDC);
357: }
358: // *
359: // * W := W * V1'
360: // *
361: this._dtrmm.Run("Right", "Lower", "Transpose", "Unit", N, K
362: , ONE, V, offset_v, LDV, ref WORK, offset_work, LDWORK);
363: // *
364: // * C1 := C1 - W'
365: // *
366: for (J = 1; J <= K; J++)
367: {
368: for (I = 1; I <= N; I++)
369: {
370: C[J+I * LDC + o_c] = C[J+I * LDC + o_c] - WORK[I+J * LDWORK + o_work];
371: }
372: }
373: // *
374: }
375: else
376: {
377: if (this._lsame.Run(SIDE, "R"))
378: {
379: // *
380: // * Form C * H or C * H' where C = ( C1 C2 )
381: // *
382: // * W := C * V = (C1*V1 + C2*V2) (stored in WORK)
383: // *
384: // * W := C1
385: // *
386: for (J = 1; J <= K; J++)
387: {
388: this._dcopy.Run(M, C, 1+J * LDC + o_c, 1, ref WORK, 1+J * LDWORK + o_work, 1);
389: }
390: // *
391: // * W := W * V1
392: // *
393: this._dtrmm.Run("Right", "Lower", "No transpose", "Unit", M, K
394: , ONE, V, offset_v, LDV, ref WORK, offset_work, LDWORK);
395: if (N > K)
396: {
397: // *
398: // * W := W + C2 * V2
399: // *
400: this._dgemm.Run("No transpose", "No transpose", M, K, N - K, ONE
401: , C, 1+(K + 1) * LDC + o_c, LDC, V, K + 1+1 * LDV + o_v, LDV, ONE, ref WORK, offset_work
402: , LDWORK);
403: }
404: // *
405: // * W := W * T or W * T'
406: // *
407: this._dtrmm.Run("Right", "Upper", TRANS, "Non-unit", M, K
408: , ONE, T, offset_t, LDT, ref WORK, offset_work, LDWORK);
409: // *
410: // * C := C - W * V'
411: // *
412: if (N > K)
413: {
414: // *
415: // * C2 := C2 - W * V2'
416: // *
417: this._dgemm.Run("No transpose", "Transpose", M, N - K, K, - ONE
418: , WORK, offset_work, LDWORK, V, K + 1+1 * LDV + o_v, LDV, ONE, ref C, 1+(K + 1) * LDC + o_c
419: , LDC);
420: }
421: // *
422: // * W := W * V1'
423: // *
424: this._dtrmm.Run("Right", "Lower", "Transpose", "Unit", M, K
425: , ONE, V, offset_v, LDV, ref WORK, offset_work, LDWORK);
426: // *
427: // * C1 := C1 - W
428: // *
429: for (J = 1; J <= K; J++)
430: {
431: for (I = 1; I <= M; I++)
432: {
433: C[I+J * LDC + o_c] = C[I+J * LDC + o_c] - WORK[I+J * LDWORK + o_work];
434: }
435: }
436: }
437: }
438: // *
439: }
440: else
441: {
442: // *
443: // * Let V = ( V1 )
444: // * ( V2 ) (last K rows)
445: // * where V2 is unit upper triangular.
446: // *
447: if (this._lsame.Run(SIDE, "L"))
448: {
449: // *
450: // * Form H * C or H' * C where C = ( C1 )
451: // * ( C2 )
452: // *
453: // * W := C' * V = (C1'*V1 + C2'*V2) (stored in WORK)
454: // *
455: // * W := C2'
456: // *
457: for (J = 1; J <= K; J++)
458: {
459: this._dcopy.Run(N, C, M - K + J+1 * LDC + o_c, LDC, ref WORK, 1+J * LDWORK + o_work, 1);
460: }
461: // *
462: // * W := W * V2
463: // *
464: this._dtrmm.Run("Right", "Upper", "No transpose", "Unit", N, K
465: , ONE, V, M - K + 1+1 * LDV + o_v, LDV, ref WORK, offset_work, LDWORK);
466: if (M > K)
467: {
468: // *
469: // * W := W + C1'*V1
470: // *
471: this._dgemm.Run("Transpose", "No transpose", N, K, M - K, ONE
472: , C, offset_c, LDC, V, offset_v, LDV, ONE, ref WORK, offset_work
473: , LDWORK);
474: }
475: // *
476: // * W := W * T' or W * T
477: // *
478: this._dtrmm.Run("Right", "Lower", TRANST, "Non-unit", N, K
479: , ONE, T, offset_t, LDT, ref WORK, offset_work, LDWORK);
480: // *
481: // * C := C - V * W'
482: // *
483: if (M > K)
484: {
485: // *
486: // * C1 := C1 - V1 * W'
487: // *
488: this._dgemm.Run("No transpose", "Transpose", M - K, N, K, - ONE
489: , V, offset_v, LDV, WORK, offset_work, LDWORK, ONE, ref C, offset_c
490: , LDC);
491: }
492: // *
493: // * W := W * V2'
494: // *
495: this._dtrmm.Run("Right", "Upper", "Transpose", "Unit", N, K
496: , ONE, V, M - K + 1+1 * LDV + o_v, LDV, ref WORK, offset_work, LDWORK);
497: // *
498: // * C2 := C2 - W'
499: // *
500: for (J = 1; J <= K; J++)
501: {
502: for (I = 1; I <= N; I++)
503: {
504: C[M - K + J+I * LDC + o_c] = C[M - K + J+I * LDC + o_c] - WORK[I+J * LDWORK + o_work];
505: }
506: }
507: // *
508: }
509: else
510: {
511: if (this._lsame.Run(SIDE, "R"))
512: {
513: // *
514: // * Form C * H or C * H' where C = ( C1 C2 )
515: // *
516: // * W := C * V = (C1*V1 + C2*V2) (stored in WORK)
517: // *
518: // * W := C2
519: // *
520: for (J = 1; J <= K; J++)
521: {
522: this._dcopy.Run(M, C, 1+(N - K + J) * LDC + o_c, 1, ref WORK, 1+J * LDWORK + o_work, 1);
523: }
524: // *
525: // * W := W * V2
526: // *
527: this._dtrmm.Run("Right", "Upper", "No transpose", "Unit", M, K
528: , ONE, V, N - K + 1+1 * LDV + o_v, LDV, ref WORK, offset_work, LDWORK);
529: if (N > K)
530: {
531: // *
532: // * W := W + C1 * V1
533: // *
534: this._dgemm.Run("No transpose", "No transpose", M, K, N - K, ONE
535: , C, offset_c, LDC, V, offset_v, LDV, ONE, ref WORK, offset_work
536: , LDWORK);
537: }
538: // *
539: // * W := W * T or W * T'
540: // *
541: this._dtrmm.Run("Right", "Lower", TRANS, "Non-unit", M, K
542: , ONE, T, offset_t, LDT, ref WORK, offset_work, LDWORK);
543: // *
544: // * C := C - W * V'
545: // *
546: if (N > K)
547: {
548: // *
549: // * C1 := C1 - W * V1'
550: // *
551: this._dgemm.Run("No transpose", "Transpose", M, N - K, K, - ONE
552: , WORK, offset_work, LDWORK, V, offset_v, LDV, ONE, ref C, offset_c
553: , LDC);
554: }
555: // *
556: // * W := W * V2'
557: // *
558: this._dtrmm.Run("Right", "Upper", "Transpose", "Unit", M, K
559: , ONE, V, N - K + 1+1 * LDV + o_v, LDV, ref WORK, offset_work, LDWORK);
560: // *
561: // * C2 := C2 - W
562: // *
563: for (J = 1; J <= K; J++)
564: {
565: for (I = 1; I <= M; I++)
566: {
567: C[I+(N - K + J) * LDC + o_c] = C[I+(N - K + J) * LDC + o_c] - WORK[I+J * LDWORK + o_work];
568: }
569: }
570: }
571: }
572: }
573: // *
574: }
575: else
576: {
577: if (this._lsame.Run(STOREV, "R"))
578: {
579: // *
580: if (this._lsame.Run(DIRECT, "F"))
581: {
582: // *
583: // * Let V = ( V1 V2 ) (V1: first K columns)
584: // * where V1 is unit upper triangular.
585: // *
586: if (this._lsame.Run(SIDE, "L"))
587: {
588: // *
589: // * Form H * C or H' * C where C = ( C1 )
590: // * ( C2 )
591: // *
592: // * W := C' * V' = (C1'*V1' + C2'*V2') (stored in WORK)
593: // *
594: // * W := C1'
595: // *
596: for (J = 1; J <= K; J++)
597: {
598: this._dcopy.Run(N, C, J+1 * LDC + o_c, LDC, ref WORK, 1+J * LDWORK + o_work, 1);
599: }
600: // *
601: // * W := W * V1'
602: // *
603: this._dtrmm.Run("Right", "Upper", "Transpose", "Unit", N, K
604: , ONE, V, offset_v, LDV, ref WORK, offset_work, LDWORK);
605: if (M > K)
606: {
607: // *
608: // * W := W + C2'*V2'
609: // *
610: this._dgemm.Run("Transpose", "Transpose", N, K, M - K, ONE
611: , C, K + 1+1 * LDC + o_c, LDC, V, 1+(K + 1) * LDV + o_v, LDV, ONE, ref WORK, offset_work
612: , LDWORK);
613: }
614: // *
615: // * W := W * T' or W * T
616: // *
617: this._dtrmm.Run("Right", "Upper", TRANST, "Non-unit", N, K
618: , ONE, T, offset_t, LDT, ref WORK, offset_work, LDWORK);
619: // *
620: // * C := C - V' * W'
621: // *
622: if (M > K)
623: {
624: // *
625: // * C2 := C2 - V2' * W'
626: // *
627: this._dgemm.Run("Transpose", "Transpose", M - K, N, K, - ONE
628: , V, 1+(K + 1) * LDV + o_v, LDV, WORK, offset_work, LDWORK, ONE, ref C, K + 1+1 * LDC + o_c
629: , LDC);
630: }
631: // *
632: // * W := W * V1
633: // *
634: this._dtrmm.Run("Right", "Upper", "No transpose", "Unit", N, K
635: , ONE, V, offset_v, LDV, ref WORK, offset_work, LDWORK);
636: // *
637: // * C1 := C1 - W'
638: // *
639: for (J = 1; J <= K; J++)
640: {
641: for (I = 1; I <= N; I++)
642: {
643: C[J+I * LDC + o_c] = C[J+I * LDC + o_c] - WORK[I+J * LDWORK + o_work];
644: }
645: }
646: // *
647: }
648: else
649: {
650: if (this._lsame.Run(SIDE, "R"))
651: {
652: // *
653: // * Form C * H or C * H' where C = ( C1 C2 )
654: // *
655: // * W := C * V' = (C1*V1' + C2*V2') (stored in WORK)
656: // *
657: // * W := C1
658: // *
659: for (J = 1; J <= K; J++)
660: {
661: this._dcopy.Run(M, C, 1+J * LDC + o_c, 1, ref WORK, 1+J * LDWORK + o_work, 1);
662: }
663: // *
664: // * W := W * V1'
665: // *
666: this._dtrmm.Run("Right", "Upper", "Transpose", "Unit", M, K
667: , ONE, V, offset_v, LDV, ref WORK, offset_work, LDWORK);
668: if (N > K)
669: {
670: // *
671: // * W := W + C2 * V2'
672: // *
673: this._dgemm.Run("No transpose", "Transpose", M, K, N - K, ONE
674: , C, 1+(K + 1) * LDC + o_c, LDC, V, 1+(K + 1) * LDV + o_v, LDV, ONE, ref WORK, offset_work
675: , LDWORK);
676: }
677: // *
678: // * W := W * T or W * T'
679: // *
680: this._dtrmm.Run("Right", "Upper", TRANS, "Non-unit", M, K
681: , ONE, T, offset_t, LDT, ref WORK, offset_work, LDWORK);
682: // *
683: // * C := C - W * V
684: // *
685: if (N > K)
686: {
687: // *
688: // * C2 := C2 - W * V2
689: // *
690: this._dgemm.Run("No transpose", "No transpose", M, N - K, K, - ONE
691: , WORK, offset_work, LDWORK, V, 1+(K + 1) * LDV + o_v, LDV, ONE, ref C, 1+(K + 1) * LDC + o_c
692: , LDC);
693: }
694: // *
695: // * W := W * V1
696: // *
697: this._dtrmm.Run("Right", "Upper", "No transpose", "Unit", M, K
698: , ONE, V, offset_v, LDV, ref WORK, offset_work, LDWORK);
699: // *
700: // * C1 := C1 - W
701: // *
702: for (J = 1; J <= K; J++)
703: {
704: for (I = 1; I <= M; I++)
705: {
706: C[I+J * LDC + o_c] = C[I+J * LDC + o_c] - WORK[I+J * LDWORK + o_work];
707: }
708: }
709: // *
710: }
711: }
712: // *
713: }
714: else
715: {
716: // *
717: // * Let V = ( V1 V2 ) (V2: last K columns)
718: // * where V2 is unit lower triangular.
719: // *
720: if (this._lsame.Run(SIDE, "L"))
721: {
722: // *
723: // * Form H * C or H' * C where C = ( C1 )
724: // * ( C2 )
725: // *
726: // * W := C' * V' = (C1'*V1' + C2'*V2') (stored in WORK)
727: // *
728: // * W := C2'
729: // *
730: for (J = 1; J <= K; J++)
731: {
732: this._dcopy.Run(N, C, M - K + J+1 * LDC + o_c, LDC, ref WORK, 1+J * LDWORK + o_work, 1);
733: }
734: // *
735: // * W := W * V2'
736: // *
737: this._dtrmm.Run("Right", "Lower", "Transpose", "Unit", N, K
738: , ONE, V, 1+(M - K + 1) * LDV + o_v, LDV, ref WORK, offset_work, LDWORK);
739: if (M > K)
740: {
741: // *
742: // * W := W + C1'*V1'
743: // *
744: this._dgemm.Run("Transpose", "Transpose", N, K, M - K, ONE
745: , C, offset_c, LDC, V, offset_v, LDV, ONE, ref WORK, offset_work
746: , LDWORK);
747: }
748: // *
749: // * W := W * T' or W * T
750: // *
751: this._dtrmm.Run("Right", "Lower", TRANST, "Non-unit", N, K
752: , ONE, T, offset_t, LDT, ref WORK, offset_work, LDWORK);
753: // *
754: // * C := C - V' * W'
755: // *
756: if (M > K)
757: {
758: // *
759: // * C1 := C1 - V1' * W'
760: // *
761: this._dgemm.Run("Transpose", "Transpose", M - K, N, K, - ONE
762: , V, offset_v, LDV, WORK, offset_work, LDWORK, ONE, ref C, offset_c
763: , LDC);
764: }
765: // *
766: // * W := W * V2
767: // *
768: this._dtrmm.Run("Right", "Lower", "No transpose", "Unit", N, K
769: , ONE, V, 1+(M - K + 1) * LDV + o_v, LDV, ref WORK, offset_work, LDWORK);
770: // *
771: // * C2 := C2 - W'
772: // *
773: for (J = 1; J <= K; J++)
774: {
775: for (I = 1; I <= N; I++)
776: {
777: C[M - K + J+I * LDC + o_c] = C[M - K + J+I * LDC + o_c] - WORK[I+J * LDWORK + o_work];
778: }
779: }
780: // *
781: }
782: else
783: {
784: if (this._lsame.Run(SIDE, "R"))
785: {
786: // *
787: // * Form C * H or C * H' where C = ( C1 C2 )
788: // *
789: // * W := C * V' = (C1*V1' + C2*V2') (stored in WORK)
790: // *
791: // * W := C2
792: // *
793: for (J = 1; J <= K; J++)
794: {
795: this._dcopy.Run(M, C, 1+(N - K + J) * LDC + o_c, 1, ref WORK, 1+J * LDWORK + o_work, 1);
796: }
797: // *
798: // * W := W * V2'
799: // *
800: this._dtrmm.Run("Right", "Lower", "Transpose", "Unit", M, K
801: , ONE, V, 1+(N - K + 1) * LDV + o_v, LDV, ref WORK, offset_work, LDWORK);
802: if (N > K)
803: {
804: // *
805: // * W := W + C1 * V1'
806: // *
807: this._dgemm.Run("No transpose", "Transpose", M, K, N - K, ONE
808: , C, offset_c, LDC, V, offset_v, LDV, ONE, ref WORK, offset_work
809: , LDWORK);
810: }
811: // *
812: // * W := W * T or W * T'
813: // *
814: this._dtrmm.Run("Right", "Lower", TRANS, "Non-unit", M, K
815: , ONE, T, offset_t, LDT, ref WORK, offset_work, LDWORK);
816: // *
817: // * C := C - W * V
818: // *
819: if (N > K)
820: {
821: // *
822: // * C1 := C1 - W * V1
823: // *
824: this._dgemm.Run("No transpose", "No transpose", M, N - K, K, - ONE
825: , WORK, offset_work, LDWORK, V, offset_v, LDV, ONE, ref C, offset_c
826: , LDC);
827: }
828: // *
829: // * W := W * V2
830: // *
831: this._dtrmm.Run("Right", "Lower", "No transpose", "Unit", M, K
832: , ONE, V, 1+(N - K + 1) * LDV + o_v, LDV, ref WORK, offset_work, LDWORK);
833: // *
834: // * C1 := C1 - W
835: // *
836: for (J = 1; J <= K; J++)
837: {
838: for (I = 1; I <= M; I++)
839: {
840: C[I+(N - K + J) * LDC + o_c] = C[I+(N - K + J) * LDC + o_c] - WORK[I+J * LDWORK + o_work];
841: }
842: }
843: // *
844: }
845: }
846: // *
847: }
848: }
849: }
850: // *
851: return;
852: // *
853: // * End of DLARFB
854: // *
855:
856: #endregion
857:
858: }
859: }
860: }