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 routine (version 3.1) --
21: /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
22: /// November 2006
23: /// Purpose
24: /// =======
25: ///
26: /// DGGSVP computes orthogonal matrices U, V and Q such that
27: ///
28: /// N-K-L K L
29: /// U'*A*Q = K ( 0 A12 A13 ) if M-K-L .GE. 0;
30: /// L ( 0 0 A23 )
31: /// M-K-L ( 0 0 0 )
32: ///
33: /// N-K-L K L
34: /// = K ( 0 A12 A13 ) if M-K-L .LT. 0;
35: /// M-K ( 0 0 A23 )
36: ///
37: /// N-K-L K L
38: /// V'*B*Q = L ( 0 0 B13 )
39: /// P-L ( 0 0 0 )
40: ///
41: /// where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular
42: /// upper triangular; A23 is L-by-L upper triangular if M-K-L .GE. 0,
43: /// otherwise A23 is (M-K)-by-L upper trapezoidal. K+L = the effective
44: /// numerical rank of the (M+P)-by-N matrix (A',B')'. Z' denotes the
45: /// transpose of Z.
46: ///
47: /// This decomposition is the preprocessing step for computing the
48: /// Generalized Singular Value Decomposition (GSVD), see subroutine
49: /// DGGSVD.
50: ///
51: ///</summary>
52: public class DGGSVP
53: {
54:
55:
56: #region Dependencies
57:
58: LSAME _lsame; DGEQPF _dgeqpf; DGEQR2 _dgeqr2; DGERQ2 _dgerq2; DLACPY _dlacpy; DLAPMT _dlapmt; DLASET _dlaset;
59: DORG2R _dorg2r;DORM2R _dorm2r; DORMR2 _dormr2; XERBLA _xerbla;
60:
61: #endregion
62:
63:
64: #region Fields
65:
66: const double ZERO = 0.0E+0; const double ONE = 1.0E+0; bool FORWRD = false; bool WANTQ = false; bool WANTU = false;
67: bool WANTV = false;int I = 0; int J = 0;
68:
69: #endregion
70:
71: public DGGSVP(LSAME lsame, DGEQPF dgeqpf, DGEQR2 dgeqr2, DGERQ2 dgerq2, DLACPY dlacpy, DLAPMT dlapmt, DLASET dlaset, DORG2R dorg2r, DORM2R dorm2r, DORMR2 dormr2
72: , XERBLA xerbla)
73: {
74:
75:
76: #region Set Dependencies
77:
78: this._lsame = lsame; this._dgeqpf = dgeqpf; this._dgeqr2 = dgeqr2; this._dgerq2 = dgerq2; this._dlacpy = dlacpy;
79: this._dlapmt = dlapmt;this._dlaset = dlaset; this._dorg2r = dorg2r; this._dorm2r = dorm2r; this._dormr2 = dormr2;
80: this._xerbla = xerbla;
81:
82: #endregion
83:
84: }
85:
86: public DGGSVP()
87: {
88:
89:
90: #region Dependencies (Initialization)
91:
92: LSAME lsame = new LSAME();
93: XERBLA xerbla = new XERBLA();
94: DLAMC3 dlamc3 = new DLAMC3();
95: DLAPY2 dlapy2 = new DLAPY2();
96: DNRM2 dnrm2 = new DNRM2();
97: DSCAL dscal = new DSCAL();
98: DSWAP dswap = new DSWAP();
99: IDAMAX idamax = new IDAMAX();
100: DLAPMT dlapmt = new DLAPMT();
101: DGEMV dgemv = new DGEMV(lsame, xerbla);
102: DGER dger = new DGER(xerbla);
103: DLARF dlarf = new DLARF(dgemv, dger, lsame);
104: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
105: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
106: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
107: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
108: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
109: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
110: DGEQR2 dgeqr2 = new DGEQR2(dlarf, dlarfg, xerbla);
111: DORM2R dorm2r = new DORM2R(lsame, dlarf, xerbla);
112: DGEQPF dgeqpf = new DGEQPF(dgeqr2, dlarf, dlarfg, dorm2r, dswap, xerbla, idamax, dlamch, dnrm2);
113: DGERQ2 dgerq2 = new DGERQ2(dlarf, dlarfg, xerbla);
114: DLACPY dlacpy = new DLACPY(lsame);
115: DLASET dlaset = new DLASET(lsame);
116: DORG2R dorg2r = new DORG2R(dlarf, dscal, xerbla);
117: DORMR2 dormr2 = new DORMR2(lsame, dlarf, xerbla);
118:
119: #endregion
120:
121:
122: #region Set Dependencies
123:
124: this._lsame = lsame; this._dgeqpf = dgeqpf; this._dgeqr2 = dgeqr2; this._dgerq2 = dgerq2; this._dlacpy = dlacpy;
125: this._dlapmt = dlapmt;this._dlaset = dlaset; this._dorg2r = dorg2r; this._dorm2r = dorm2r; this._dormr2 = dormr2;
126: this._xerbla = xerbla;
127:
128: #endregion
129:
130: }
131: /// <summary>
132: /// Purpose
133: /// =======
134: ///
135: /// DGGSVP computes orthogonal matrices U, V and Q such that
136: ///
137: /// N-K-L K L
138: /// U'*A*Q = K ( 0 A12 A13 ) if M-K-L .GE. 0;
139: /// L ( 0 0 A23 )
140: /// M-K-L ( 0 0 0 )
141: ///
142: /// N-K-L K L
143: /// = K ( 0 A12 A13 ) if M-K-L .LT. 0;
144: /// M-K ( 0 0 A23 )
145: ///
146: /// N-K-L K L
147: /// V'*B*Q = L ( 0 0 B13 )
148: /// P-L ( 0 0 0 )
149: ///
150: /// where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular
151: /// upper triangular; A23 is L-by-L upper triangular if M-K-L .GE. 0,
152: /// otherwise A23 is (M-K)-by-L upper trapezoidal. K+L = the effective
153: /// numerical rank of the (M+P)-by-N matrix (A',B')'. Z' denotes the
154: /// transpose of Z.
155: ///
156: /// This decomposition is the preprocessing step for computing the
157: /// Generalized Singular Value Decomposition (GSVD), see subroutine
158: /// DGGSVD.
159: ///
160: ///</summary>
161: /// <param name="JOBU">
162: /// (input) CHARACTER*1
163: /// = 'U': Orthogonal matrix U is computed;
164: /// = 'N': U is not computed.
165: ///</param>
166: /// <param name="JOBV">
167: /// (input) CHARACTER*1
168: /// = 'V': Orthogonal matrix V is computed;
169: /// = 'N': V is not computed.
170: ///</param>
171: /// <param name="JOBQ">
172: /// (input) CHARACTER*1
173: /// = 'Q': Orthogonal matrix Q is computed;
174: /// = 'N': Q is not computed.
175: ///</param>
176: /// <param name="M">
177: /// (input) INTEGER
178: /// The number of rows of the matrix A. M .GE. 0.
179: ///</param>
180: /// <param name="P">
181: /// (input) INTEGER
182: /// The number of rows of the matrix B. P .GE. 0.
183: ///</param>
184: /// <param name="N">
185: /// (input) INTEGER
186: /// The number of columns of the matrices A and B. N .GE. 0.
187: ///</param>
188: /// <param name="A">
189: /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
190: /// On entry, the M-by-N matrix A.
191: /// On exit, A contains the triangular (or trapezoidal) matrix
192: /// described in the Purpose section.
193: ///</param>
194: /// <param name="LDA">
195: /// (input) INTEGER
196: /// The leading dimension of the array A. LDA .GE. max(1,M).
197: ///</param>
198: /// <param name="B">
199: /// (input/output) DOUBLE PRECISION array, dimension (LDB,N)
200: /// On entry, the P-by-N matrix B.
201: /// On exit, B contains the triangular matrix described in
202: /// the Purpose section.
203: ///</param>
204: /// <param name="LDB">
205: /// (input) INTEGER
206: /// The leading dimension of the array B. LDB .GE. max(1,P).
207: ///</param>
208: /// <param name="TOLA">
209: /// (input) DOUBLE PRECISION
210: ///</param>
211: /// <param name="TOLB">
212: /// (input) DOUBLE PRECISION
213: /// TOLA and TOLB are the thresholds to determine the effective
214: /// numerical rank of matrix B and a subblock of A. Generally,
215: /// they are set to
216: /// TOLA = MAX(M,N)*norm(A)*MAZHEPS,
217: /// TOLB = MAX(P,N)*norm(B)*MAZHEPS.
218: /// The size of TOLA and TOLB may affect the size of backward
219: /// errors of the decomposition.
220: ///</param>
221: /// <param name="K">
222: /// (output) INTEGER
223: ///</param>
224: /// <param name="L">
225: /// ( 0 0 A23 )
226: ///</param>
227: /// <param name="U">
228: /// (output) DOUBLE PRECISION array, dimension (LDU,M)
229: /// If JOBU = 'U', U contains the orthogonal matrix U.
230: /// If JOBU = 'N', U is not referenced.
231: ///</param>
232: /// <param name="LDU">
233: /// (input) INTEGER
234: /// The leading dimension of the array U. LDU .GE. max(1,M) if
235: /// JOBU = 'U'; LDU .GE. 1 otherwise.
236: ///</param>
237: /// <param name="V">
238: /// (output) DOUBLE PRECISION array, dimension (LDV,M)
239: /// If JOBV = 'V', V contains the orthogonal matrix V.
240: /// If JOBV = 'N', V is not referenced.
241: ///</param>
242: /// <param name="LDV">
243: /// (input) INTEGER
244: /// The leading dimension of the array V. LDV .GE. max(1,P) if
245: /// JOBV = 'V'; LDV .GE. 1 otherwise.
246: ///</param>
247: /// <param name="Q">
248: /// (output) DOUBLE PRECISION array, dimension (LDQ,N)
249: /// If JOBQ = 'Q', Q contains the orthogonal matrix Q.
250: /// If JOBQ = 'N', Q is not referenced.
251: ///</param>
252: /// <param name="LDQ">
253: /// (input) INTEGER
254: /// The leading dimension of the array Q. LDQ .GE. max(1,N) if
255: /// JOBQ = 'Q'; LDQ .GE. 1 otherwise.
256: ///</param>
257: /// <param name="IWORK">
258: /// (workspace) INTEGER array, dimension (N)
259: ///</param>
260: /// <param name="TAU">
261: /// (workspace) DOUBLE PRECISION array, dimension (N)
262: ///</param>
263: /// <param name="WORK">
264: /// (workspace) DOUBLE PRECISION array, dimension (max(3*N,M,P))
265: ///</param>
266: /// <param name="INFO">
267: /// (output) INTEGER
268: /// = 0: successful exit
269: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
270: ///
271: ///</param>
272: public void Run(string JOBU, string JOBV, string JOBQ, int M, int P, int N
273: , ref double[] A, int offset_a, int LDA, ref double[] B, int offset_b, int LDB, double TOLA, double TOLB
274: , ref int K, ref int L, ref double[] U, int offset_u, int LDU, ref double[] V, int offset_v, int LDV
275: , ref double[] Q, int offset_q, int LDQ, ref int[] IWORK, int offset_iwork, ref double[] TAU, int offset_tau, ref double[] WORK, int offset_work, ref int INFO)
276: {
277:
278: #region Array Index Correction
279:
280: int o_a = -1 - LDA + offset_a; int o_b = -1 - LDB + offset_b; int o_u = -1 - LDU + offset_u;
281: int o_v = -1 - LDV + offset_v; int o_q = -1 - LDQ + offset_q; int o_iwork = -1 + offset_iwork;
282: int o_tau = -1 + offset_tau; int o_work = -1 + offset_work;
283:
284: #endregion
285:
286:
287: #region Strings
288:
289: JOBU = JOBU.Substring(0, 1); JOBV = JOBV.Substring(0, 1); JOBQ = JOBQ.Substring(0, 1);
290:
291: #endregion
292:
293:
294: #region Prolog
295:
296: // *
297: // * -- LAPACK routine (version 3.1) --
298: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
299: // * November 2006
300: // *
301: // * .. Scalar Arguments ..
302: // * ..
303: // * .. Array Arguments ..
304: // * ..
305: // *
306: // * Purpose
307: // * =======
308: // *
309: // * DGGSVP computes orthogonal matrices U, V and Q such that
310: // *
311: // * N-K-L K L
312: // * U'*A*Q = K ( 0 A12 A13 ) if M-K-L >= 0;
313: // * L ( 0 0 A23 )
314: // * M-K-L ( 0 0 0 )
315: // *
316: // * N-K-L K L
317: // * = K ( 0 A12 A13 ) if M-K-L < 0;
318: // * M-K ( 0 0 A23 )
319: // *
320: // * N-K-L K L
321: // * V'*B*Q = L ( 0 0 B13 )
322: // * P-L ( 0 0 0 )
323: // *
324: // * where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular
325: // * upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0,
326: // * otherwise A23 is (M-K)-by-L upper trapezoidal. K+L = the effective
327: // * numerical rank of the (M+P)-by-N matrix (A',B')'. Z' denotes the
328: // * transpose of Z.
329: // *
330: // * This decomposition is the preprocessing step for computing the
331: // * Generalized Singular Value Decomposition (GSVD), see subroutine
332: // * DGGSVD.
333: // *
334: // * Arguments
335: // * =========
336: // *
337: // * JOBU (input) CHARACTER*1
338: // * = 'U': Orthogonal matrix U is computed;
339: // * = 'N': U is not computed.
340: // *
341: // * JOBV (input) CHARACTER*1
342: // * = 'V': Orthogonal matrix V is computed;
343: // * = 'N': V is not computed.
344: // *
345: // * JOBQ (input) CHARACTER*1
346: // * = 'Q': Orthogonal matrix Q is computed;
347: // * = 'N': Q is not computed.
348: // *
349: // * M (input) INTEGER
350: // * The number of rows of the matrix A. M >= 0.
351: // *
352: // * P (input) INTEGER
353: // * The number of rows of the matrix B. P >= 0.
354: // *
355: // * N (input) INTEGER
356: // * The number of columns of the matrices A and B. N >= 0.
357: // *
358: // * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
359: // * On entry, the M-by-N matrix A.
360: // * On exit, A contains the triangular (or trapezoidal) matrix
361: // * described in the Purpose section.
362: // *
363: // * LDA (input) INTEGER
364: // * The leading dimension of the array A. LDA >= max(1,M).
365: // *
366: // * B (input/output) DOUBLE PRECISION array, dimension (LDB,N)
367: // * On entry, the P-by-N matrix B.
368: // * On exit, B contains the triangular matrix described in
369: // * the Purpose section.
370: // *
371: // * LDB (input) INTEGER
372: // * The leading dimension of the array B. LDB >= max(1,P).
373: // *
374: // * TOLA (input) DOUBLE PRECISION
375: // * TOLB (input) DOUBLE PRECISION
376: // * TOLA and TOLB are the thresholds to determine the effective
377: // * numerical rank of matrix B and a subblock of A. Generally,
378: // * they are set to
379: // * TOLA = MAX(M,N)*norm(A)*MAZHEPS,
380: // * TOLB = MAX(P,N)*norm(B)*MAZHEPS.
381: // * The size of TOLA and TOLB may affect the size of backward
382: // * errors of the decomposition.
383: // *
384: // * K (output) INTEGER
385: // * L (output) INTEGER
386: // * On exit, K and L specify the dimension of the subblocks
387: // * described in Purpose.
388: // * K + L = effective numerical rank of (A',B')'.
389: // *
390: // * U (output) DOUBLE PRECISION array, dimension (LDU,M)
391: // * If JOBU = 'U', U contains the orthogonal matrix U.
392: // * If JOBU = 'N', U is not referenced.
393: // *
394: // * LDU (input) INTEGER
395: // * The leading dimension of the array U. LDU >= max(1,M) if
396: // * JOBU = 'U'; LDU >= 1 otherwise.
397: // *
398: // * V (output) DOUBLE PRECISION array, dimension (LDV,M)
399: // * If JOBV = 'V', V contains the orthogonal matrix V.
400: // * If JOBV = 'N', V is not referenced.
401: // *
402: // * LDV (input) INTEGER
403: // * The leading dimension of the array V. LDV >= max(1,P) if
404: // * JOBV = 'V'; LDV >= 1 otherwise.
405: // *
406: // * Q (output) DOUBLE PRECISION array, dimension (LDQ,N)
407: // * If JOBQ = 'Q', Q contains the orthogonal matrix Q.
408: // * If JOBQ = 'N', Q is not referenced.
409: // *
410: // * LDQ (input) INTEGER
411: // * The leading dimension of the array Q. LDQ >= max(1,N) if
412: // * JOBQ = 'Q'; LDQ >= 1 otherwise.
413: // *
414: // * IWORK (workspace) INTEGER array, dimension (N)
415: // *
416: // * TAU (workspace) DOUBLE PRECISION array, dimension (N)
417: // *
418: // * WORK (workspace) DOUBLE PRECISION array, dimension (max(3*N,M,P))
419: // *
420: // * INFO (output) INTEGER
421: // * = 0: successful exit
422: // * < 0: if INFO = -i, the i-th argument had an illegal value.
423: // *
424: // *
425: // * Further Details
426: // * ===============
427: // *
428: // * The subroutine uses LAPACK subroutine DGEQPF for the QR factorization
429: // * with column pivoting to detect the effective numerical rank of the
430: // * a matrix. It may be replaced by a better rank determination strategy.
431: // *
432: // * =====================================================================
433: // *
434: // * .. Parameters ..
435: // * ..
436: // * .. Local Scalars ..
437: // * ..
438: // * .. External Functions ..
439: // * ..
440: // * .. External Subroutines ..
441: // * ..
442: // * .. Intrinsic Functions ..
443: // INTRINSIC ABS, MAX, MIN;
444: // * ..
445: // * .. Executable Statements ..
446: // *
447: // * Test the input parameters
448: // *
449:
450: #endregion
451:
452:
453: #region Body
454:
455: WANTU = this._lsame.Run(JOBU, "U");
456: WANTV = this._lsame.Run(JOBV, "V");
457: WANTQ = this._lsame.Run(JOBQ, "Q");
458: FORWRD = true;
459: // *
460: INFO = 0;
461: if (!(WANTU || this._lsame.Run(JOBU, "N")))
462: {
463: INFO = - 1;
464: }
465: else
466: {
467: if (!(WANTV || this._lsame.Run(JOBV, "N")))
468: {
469: INFO = - 2;
470: }
471: else
472: {
473: if (!(WANTQ || this._lsame.Run(JOBQ, "N")))
474: {
475: INFO = - 3;
476: }
477: else
478: {
479: if (M < 0)
480: {
481: INFO = - 4;
482: }
483: else
484: {
485: if (P < 0)
486: {
487: INFO = - 5;
488: }
489: else
490: {
491: if (N < 0)
492: {
493: INFO = - 6;
494: }
495: else
496: {
497: if (LDA < Math.Max(1, M))
498: {
499: INFO = - 8;
500: }
501: else
502: {
503: if (LDB < Math.Max(1, P))
504: {
505: INFO = - 10;
506: }
507: else
508: {
509: if (LDU < 1 || (WANTU && LDU < M))
510: {
511: INFO = - 16;
512: }
513: else
514: {
515: if (LDV < 1 || (WANTV && LDV < P))
516: {
517: INFO = - 18;
518: }
519: else
520: {
521: if (LDQ < 1 || (WANTQ && LDQ < N))
522: {
523: INFO = - 20;
524: }
525: }
526: }
527: }
528: }
529: }
530: }
531: }
532: }
533: }
534: }
535: if (INFO != 0)
536: {
537: this._xerbla.Run("DGGSVP", - INFO);
538: return;
539: }
540: // *
541: // * QR with column pivoting of B: B*P = V*( S11 S12 )
542: // * ( 0 0 )
543: // *
544: for (I = 1; I <= N; I++)
545: {
546: IWORK[I + o_iwork] = 0;
547: }
548: this._dgeqpf.Run(P, N, ref B, offset_b, LDB, ref IWORK, offset_iwork, ref TAU, offset_tau
549: , ref WORK, offset_work, ref INFO);
550: // *
551: // * Update A := A*P
552: // *
553: this._dlapmt.Run(FORWRD, M, N, ref A, offset_a, LDA, ref IWORK, offset_iwork);
554: // *
555: // * Determine the effective rank of matrix B.
556: // *
557: L = 0;
558: for (I = 1; I <= Math.Min(P, N); I++)
559: {
560: if (Math.Abs(B[I+I * LDB + o_b]) > TOLB) L = L + 1;
561: }
562: // *
563: if (WANTV)
564: {
565: // *
566: // * Copy the details of V, and form V.
567: // *
568: this._dlaset.Run("Full", P, P, ZERO, ZERO, ref V, offset_v
569: , LDV);
570: if (P > 1)
571: {
572: this._dlacpy.Run("Lower", P - 1, N, B, 2+1 * LDB + o_b, LDB, ref V, 2+1 * LDV + o_v
573: , LDV);
574: }
575: this._dorg2r.Run(P, P, Math.Min(P, N), ref V, offset_v, LDV, TAU, offset_tau
576: , ref WORK, offset_work, ref INFO);
577: }
578: // *
579: // * Clean up B
580: // *
581: for (J = 1; J <= L - 1; J++)
582: {
583: for (I = J + 1; I <= L; I++)
584: {
585: B[I+J * LDB + o_b] = ZERO;
586: }
587: }
588: if (P > L)
589: {
590: this._dlaset.Run("Full", P - L, N, ZERO, ZERO, ref B, L + 1+1 * LDB + o_b
591: , LDB);
592: }
593: // *
594: if (WANTQ)
595: {
596: // *
597: // * Set Q = I and Update Q := Q*P
598: // *
599: this._dlaset.Run("Full", N, N, ZERO, ONE, ref Q, offset_q
600: , LDQ);
601: this._dlapmt.Run(FORWRD, N, N, ref Q, offset_q, LDQ, ref IWORK, offset_iwork);
602: }
603: // *
604: if (P >= L && N != L)
605: {
606: // *
607: // * RQ factorization of (S11 S12): ( S11 S12 ) = ( 0 S12 )*Z
608: // *
609: this._dgerq2.Run(L, N, ref B, offset_b, LDB, ref TAU, offset_tau, ref WORK, offset_work
610: , ref INFO);
611: // *
612: // * Update A := A*Z'
613: // *
614: this._dormr2.Run("Right", "Transpose", M, N, L, ref B, offset_b
615: , LDB, TAU, offset_tau, ref A, offset_a, LDA, ref WORK, offset_work, ref INFO);
616: // *
617: if (WANTQ)
618: {
619: // *
620: // * Update Q := Q*Z'
621: // *
622: this._dormr2.Run("Right", "Transpose", N, N, L, ref B, offset_b
623: , LDB, TAU, offset_tau, ref Q, offset_q, LDQ, ref WORK, offset_work, ref INFO);
624: }
625: // *
626: // * Clean up B
627: // *
628: this._dlaset.Run("Full", L, N - L, ZERO, ZERO, ref B, offset_b
629: , LDB);
630: for (J = N - L + 1; J <= N; J++)
631: {
632: for (I = J - N + L + 1; I <= L; I++)
633: {
634: B[I+J * LDB + o_b] = ZERO;
635: }
636: }
637: // *
638: }
639: // *
640: // * Let N-L L
641: // * A = ( A11 A12 ) M,
642: // *
643: // * then the following does the complete QR decomposition of A11:
644: // *
645: // * A11 = U*( 0 T12 )*P1'
646: // * ( 0 0 )
647: // *
648: for (I = 1; I <= N - L; I++)
649: {
650: IWORK[I + o_iwork] = 0;
651: }
652: this._dgeqpf.Run(M, N - L, ref A, offset_a, LDA, ref IWORK, offset_iwork, ref TAU, offset_tau
653: , ref WORK, offset_work, ref INFO);
654: // *
655: // * Determine the effective rank of A11
656: // *
657: K = 0;
658: for (I = 1; I <= Math.Min(M, N - L); I++)
659: {
660: if (Math.Abs(A[I+I * LDA + o_a]) > TOLA) K = K + 1;
661: }
662: // *
663: // * Update A12 := U'*A12, where A12 = A( 1:M, N-L+1:N )
664: // *
665: this._dorm2r.Run("Left", "Transpose", M, L, Math.Min(M, N - L), ref A, offset_a
666: , LDA, TAU, offset_tau, ref A, 1+(N - L + 1) * LDA + o_a, LDA, ref WORK, offset_work, ref INFO);
667: // *
668: if (WANTU)
669: {
670: // *
671: // * Copy the details of U, and form U
672: // *
673: this._dlaset.Run("Full", M, M, ZERO, ZERO, ref U, offset_u
674: , LDU);
675: if (M > 1)
676: {
677: this._dlacpy.Run("Lower", M - 1, N - L, A, 2+1 * LDA + o_a, LDA, ref U, 2+1 * LDU + o_u
678: , LDU);
679: }
680: this._dorg2r.Run(M, M, Math.Min(M, N - L), ref U, offset_u, LDU, TAU, offset_tau
681: , ref WORK, offset_work, ref INFO);
682: }
683: // *
684: if (WANTQ)
685: {
686: // *
687: // * Update Q( 1:N, 1:N-L ) = Q( 1:N, 1:N-L )*P1
688: // *
689: this._dlapmt.Run(FORWRD, N, N - L, ref Q, offset_q, LDQ, ref IWORK, offset_iwork);
690: }
691: // *
692: // * Clean up A: set the strictly lower triangular part of
693: // * A(1:K, 1:K) = 0, and A( K+1:M, 1:N-L ) = 0.
694: // *
695: for (J = 1; J <= K - 1; J++)
696: {
697: for (I = J + 1; I <= K; I++)
698: {
699: A[I+J * LDA + o_a] = ZERO;
700: }
701: }
702: if (M > K)
703: {
704: this._dlaset.Run("Full", M - K, N - L, ZERO, ZERO, ref A, K + 1+1 * LDA + o_a
705: , LDA);
706: }
707: // *
708: if (N - L > K)
709: {
710: // *
711: // * RQ factorization of ( T11 T12 ) = ( 0 T12 )*Z1
712: // *
713: this._dgerq2.Run(K, N - L, ref A, offset_a, LDA, ref TAU, offset_tau, ref WORK, offset_work
714: , ref INFO);
715: // *
716: if (WANTQ)
717: {
718: // *
719: // * Update Q( 1:N,1:N-L ) = Q( 1:N,1:N-L )*Z1'
720: // *
721: this._dormr2.Run("Right", "Transpose", N, N - L, K, ref A, offset_a
722: , LDA, TAU, offset_tau, ref Q, offset_q, LDQ, ref WORK, offset_work, ref INFO);
723: }
724: // *
725: // * Clean up A
726: // *
727: this._dlaset.Run("Full", K, N - L - K, ZERO, ZERO, ref A, offset_a
728: , LDA);
729: for (J = N - L - K + 1; J <= N - L; J++)
730: {
731: for (I = J - N + L + K + 1; I <= K; I++)
732: {
733: A[I+J * LDA + o_a] = ZERO;
734: }
735: }
736: // *
737: }
738: // *
739: if (M > K)
740: {
741: // *
742: // * QR factorization of A( K+1:M,N-L+1:N )
743: // *
744: this._dgeqr2.Run(M - K, L, ref A, K + 1+(N - L + 1) * LDA + o_a, LDA, ref TAU, offset_tau, ref WORK, offset_work
745: , ref INFO);
746: // *
747: if (WANTU)
748: {
749: // *
750: // * Update U(:,K+1:M) := U(:,K+1:M)*U1
751: // *
752: this._dorm2r.Run("Right", "No transpose", M, M - K, Math.Min(M - K, L), ref A, K + 1+(N - L + 1) * LDA + o_a
753: , LDA, TAU, offset_tau, ref U, 1+(K + 1) * LDU + o_u, LDU, ref WORK, offset_work, ref INFO);
754: }
755: // *
756: // * Clean up
757: // *
758: for (J = N - L + 1; J <= N; J++)
759: {
760: for (I = J - N + K + L + 1; I <= M; I++)
761: {
762: A[I+J * LDA + o_a] = ZERO;
763: }
764: }
765: // *
766: }
767: // *
768: return;
769: // *
770: // * End of DGGSVP
771: // *
772:
773: #endregion
774:
775: }
776: }
777: }