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: /// DTGSJA computes the generalized singular value decomposition (GSVD)
27: /// of two real upper triangular (or trapezoidal) matrices A and B.
28: ///
29: /// On entry, it is assumed that matrices A and B have the following
30: /// forms, which may be obtained by the preprocessing subroutine DGGSVP
31: /// from a general M-by-N matrix A and P-by-N matrix B:
32: ///
33: /// N-K-L K L
34: /// A = K ( 0 A12 A13 ) if M-K-L .GE. 0;
35: /// L ( 0 0 A23 )
36: /// M-K-L ( 0 0 0 )
37: ///
38: /// N-K-L K L
39: /// A = K ( 0 A12 A13 ) if M-K-L .LT. 0;
40: /// M-K ( 0 0 A23 )
41: ///
42: /// N-K-L K L
43: /// B = L ( 0 0 B13 )
44: /// P-L ( 0 0 0 )
45: ///
46: /// where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular
47: /// upper triangular; A23 is L-by-L upper triangular if M-K-L .GE. 0,
48: /// otherwise A23 is (M-K)-by-L upper trapezoidal.
49: ///
50: /// On exit,
51: ///
52: /// U'*A*Q = D1*( 0 R ), V'*B*Q = D2*( 0 R ),
53: ///
54: /// where U, V and Q are orthogonal matrices, Z' denotes the transpose
55: /// of Z, R is a nonsingular upper triangular matrix, and D1 and D2 are
56: /// ``diagonal'' matrices, which are of the following structures:
57: ///
58: /// If M-K-L .GE. 0,
59: ///
60: /// K L
61: /// D1 = K ( I 0 )
62: /// L ( 0 C )
63: /// M-K-L ( 0 0 )
64: ///
65: /// K L
66: /// D2 = L ( 0 S )
67: /// P-L ( 0 0 )
68: ///
69: /// N-K-L K L
70: /// ( 0 R ) = K ( 0 R11 R12 ) K
71: /// L ( 0 0 R22 ) L
72: ///
73: /// where
74: ///
75: /// C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
76: /// S = diag( BETA(K+1), ... , BETA(K+L) ),
77: /// C**2 + S**2 = I.
78: ///
79: /// R is stored in A(1:K+L,N-K-L+1:N) on exit.
80: ///
81: /// If M-K-L .LT. 0,
82: ///
83: /// K M-K K+L-M
84: /// D1 = K ( I 0 0 )
85: /// M-K ( 0 C 0 )
86: ///
87: /// K M-K K+L-M
88: /// D2 = M-K ( 0 S 0 )
89: /// K+L-M ( 0 0 I )
90: /// P-L ( 0 0 0 )
91: ///
92: /// N-K-L K M-K K+L-M
93: /// ( 0 R ) = K ( 0 R11 R12 R13 )
94: /// M-K ( 0 0 R22 R23 )
95: /// K+L-M ( 0 0 0 R33 )
96: ///
97: /// where
98: /// C = diag( ALPHA(K+1), ... , ALPHA(M) ),
99: /// S = diag( BETA(K+1), ... , BETA(M) ),
100: /// C**2 + S**2 = I.
101: ///
102: /// R = ( R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N) and R33 is stored
103: /// ( 0 R22 R23 )
104: /// in B(M-K+1:L,N+M-K-L+1:N) on exit.
105: ///
106: /// The computation of the orthogonal transformation matrices U, V or Q
107: /// is optional. These matrices may either be formed explicitly, or they
108: /// may be postmultiplied into input matrices U1, V1, or Q1.
109: ///
110: ///</summary>
111: public class DTGSJA
112: {
113:
114:
115: #region Dependencies
116:
117: LSAME _lsame; DCOPY _dcopy; DLAGS2 _dlags2; DLAPLL _dlapll; DLARTG _dlartg; DLASET _dlaset; DROT _drot; DSCAL _dscal;
118: XERBLA _xerbla;
119:
120: #endregion
121:
122:
123: #region Fields
124:
125: const int MAXIT = 40; const double ZERO = 0.0E+0; const double ONE = 1.0E+0; bool INITQ = false; bool INITU = false;
126: bool INITV = false;bool UPPER = false; bool WANTQ = false; bool WANTU = false; bool WANTV = false; int I = 0; int J = 0;
127: int KCYCLE = 0;double A1 = 0; double A2 = 0; double A3 = 0; double B1 = 0; double B2 = 0; double B3 = 0; double CSQ = 0;
128: double CSU = 0;double CSV = 0; double ERROR = 0; double GAMMA = 0; double RWK = 0; double SNQ = 0; double SNU = 0;
129: double SNV = 0;double SSMIN = 0;
130:
131: #endregion
132:
133: public DTGSJA(LSAME lsame, DCOPY dcopy, DLAGS2 dlags2, DLAPLL dlapll, DLARTG dlartg, DLASET dlaset, DROT drot, DSCAL dscal, XERBLA xerbla)
134: {
135:
136:
137: #region Set Dependencies
138:
139: this._lsame = lsame; this._dcopy = dcopy; this._dlags2 = dlags2; this._dlapll = dlapll; this._dlartg = dlartg;
140: this._dlaset = dlaset;this._drot = drot; this._dscal = dscal; this._xerbla = xerbla;
141:
142: #endregion
143:
144: }
145:
146: public DTGSJA()
147: {
148:
149:
150: #region Dependencies (Initialization)
151:
152: LSAME lsame = new LSAME();
153: DCOPY dcopy = new DCOPY();
154: DLAMC3 dlamc3 = new DLAMC3();
155: DDOT ddot = new DDOT();
156: DAXPY daxpy = new DAXPY();
157: DLAPY2 dlapy2 = new DLAPY2();
158: DNRM2 dnrm2 = new DNRM2();
159: DSCAL dscal = new DSCAL();
160: DLAS2 dlas2 = new DLAS2();
161: DROT drot = new DROT();
162: XERBLA xerbla = new XERBLA();
163: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
164: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
165: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
166: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
167: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
168: DLARTG dlartg = new DLARTG(dlamch);
169: DLASV2 dlasv2 = new DLASV2(dlamch);
170: DLAGS2 dlags2 = new DLAGS2(dlartg, dlasv2);
171: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
172: DLAPLL dlapll = new DLAPLL(ddot, daxpy, dlarfg, dlas2);
173: DLASET dlaset = new DLASET(lsame);
174:
175: #endregion
176:
177:
178: #region Set Dependencies
179:
180: this._lsame = lsame; this._dcopy = dcopy; this._dlags2 = dlags2; this._dlapll = dlapll; this._dlartg = dlartg;
181: this._dlaset = dlaset;this._drot = drot; this._dscal = dscal; this._xerbla = xerbla;
182:
183: #endregion
184:
185: }
186: /// <summary>
187: /// Purpose
188: /// =======
189: ///
190: /// DTGSJA computes the generalized singular value decomposition (GSVD)
191: /// of two real upper triangular (or trapezoidal) matrices A and B.
192: ///
193: /// On entry, it is assumed that matrices A and B have the following
194: /// forms, which may be obtained by the preprocessing subroutine DGGSVP
195: /// from a general M-by-N matrix A and P-by-N matrix B:
196: ///
197: /// N-K-L K L
198: /// A = K ( 0 A12 A13 ) if M-K-L .GE. 0;
199: /// L ( 0 0 A23 )
200: /// M-K-L ( 0 0 0 )
201: ///
202: /// N-K-L K L
203: /// A = K ( 0 A12 A13 ) if M-K-L .LT. 0;
204: /// M-K ( 0 0 A23 )
205: ///
206: /// N-K-L K L
207: /// B = L ( 0 0 B13 )
208: /// P-L ( 0 0 0 )
209: ///
210: /// where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular
211: /// upper triangular; A23 is L-by-L upper triangular if M-K-L .GE. 0,
212: /// otherwise A23 is (M-K)-by-L upper trapezoidal.
213: ///
214: /// On exit,
215: ///
216: /// U'*A*Q = D1*( 0 R ), V'*B*Q = D2*( 0 R ),
217: ///
218: /// where U, V and Q are orthogonal matrices, Z' denotes the transpose
219: /// of Z, R is a nonsingular upper triangular matrix, and D1 and D2 are
220: /// ``diagonal'' matrices, which are of the following structures:
221: ///
222: /// If M-K-L .GE. 0,
223: ///
224: /// K L
225: /// D1 = K ( I 0 )
226: /// L ( 0 C )
227: /// M-K-L ( 0 0 )
228: ///
229: /// K L
230: /// D2 = L ( 0 S )
231: /// P-L ( 0 0 )
232: ///
233: /// N-K-L K L
234: /// ( 0 R ) = K ( 0 R11 R12 ) K
235: /// L ( 0 0 R22 ) L
236: ///
237: /// where
238: ///
239: /// C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
240: /// S = diag( BETA(K+1), ... , BETA(K+L) ),
241: /// C**2 + S**2 = I.
242: ///
243: /// R is stored in A(1:K+L,N-K-L+1:N) on exit.
244: ///
245: /// If M-K-L .LT. 0,
246: ///
247: /// K M-K K+L-M
248: /// D1 = K ( I 0 0 )
249: /// M-K ( 0 C 0 )
250: ///
251: /// K M-K K+L-M
252: /// D2 = M-K ( 0 S 0 )
253: /// K+L-M ( 0 0 I )
254: /// P-L ( 0 0 0 )
255: ///
256: /// N-K-L K M-K K+L-M
257: /// ( 0 R ) = K ( 0 R11 R12 R13 )
258: /// M-K ( 0 0 R22 R23 )
259: /// K+L-M ( 0 0 0 R33 )
260: ///
261: /// where
262: /// C = diag( ALPHA(K+1), ... , ALPHA(M) ),
263: /// S = diag( BETA(K+1), ... , BETA(M) ),
264: /// C**2 + S**2 = I.
265: ///
266: /// R = ( R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N) and R33 is stored
267: /// ( 0 R22 R23 )
268: /// in B(M-K+1:L,N+M-K-L+1:N) on exit.
269: ///
270: /// The computation of the orthogonal transformation matrices U, V or Q
271: /// is optional. These matrices may either be formed explicitly, or they
272: /// may be postmultiplied into input matrices U1, V1, or Q1.
273: ///
274: ///</summary>
275: /// <param name="JOBU">
276: /// (input) CHARACTER*1
277: /// = 'U': U must contain an orthogonal matrix U1 on entry, and
278: /// the product U1*U is returned;
279: /// = 'I': U is initialized to the unit matrix, and the
280: /// orthogonal matrix U is returned;
281: /// = 'N': U is not computed.
282: ///</param>
283: /// <param name="JOBV">
284: /// (input) CHARACTER*1
285: /// = 'V': V must contain an orthogonal matrix V1 on entry, and
286: /// the product V1*V is returned;
287: /// = 'I': V is initialized to the unit matrix, and the
288: /// orthogonal matrix V is returned;
289: /// = 'N': V is not computed.
290: ///</param>
291: /// <param name="JOBQ">
292: /// (input) CHARACTER*1
293: /// = 'Q': Q must contain an orthogonal matrix Q1 on entry, and
294: /// the product Q1*Q is returned;
295: /// = 'I': Q is initialized to the unit matrix, and the
296: /// orthogonal matrix Q is returned;
297: /// = 'N': Q is not computed.
298: ///</param>
299: /// <param name="M">
300: /// (input) INTEGER
301: /// The number of rows of the matrix A. M .GE. 0.
302: ///</param>
303: /// <param name="P">
304: /// (input) INTEGER
305: /// The number of rows of the matrix B. P .GE. 0.
306: ///</param>
307: /// <param name="N">
308: /// (input) INTEGER
309: /// The number of columns of the matrices A and B. N .GE. 0.
310: ///</param>
311: /// <param name="K">
312: /// L
313: ///</param>
314: /// <param name="L">
315: /// ( 0 0 A23 )
316: ///</param>
317: /// <param name="A">
318: /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
319: /// On entry, the M-by-N matrix A.
320: /// On exit, A(N-K+1:N,1:MIN(K+L,M) ) contains the triangular
321: /// matrix R or part of R. See Purpose for details.
322: ///</param>
323: /// <param name="LDA">
324: /// (input) INTEGER
325: /// The leading dimension of the array A. LDA .GE. max(1,M).
326: ///</param>
327: /// <param name="B">
328: /// (input/output) DOUBLE PRECISION array, dimension (LDB,N)
329: /// On entry, the P-by-N matrix B.
330: /// On exit, if necessary, B(M-K+1:L,N+M-K-L+1:N) contains
331: /// a part of R. See Purpose for details.
332: ///</param>
333: /// <param name="LDB">
334: /// (input) INTEGER
335: /// The leading dimension of the array B. LDB .GE. max(1,P).
336: ///</param>
337: /// <param name="TOLA">
338: /// (input) DOUBLE PRECISION
339: ///</param>
340: /// <param name="TOLB">
341: /// (input) DOUBLE PRECISION
342: /// TOLA and TOLB are the convergence criteria for the Jacobi-
343: /// Kogbetliantz iteration procedure. Generally, they are the
344: /// same as used in the preprocessing step, say
345: /// TOLA = max(M,N)*norm(A)*MAZHEPS,
346: /// TOLB = max(P,N)*norm(B)*MAZHEPS.
347: ///</param>
348: /// <param name="ALPHA">
349: /// (output) DOUBLE PRECISION array, dimension (N)
350: ///</param>
351: /// <param name="BETA">
352: /// (output) DOUBLE PRECISION array, dimension (N)
353: /// On exit, ALPHA and BETA contain the generalized singular
354: /// value pairs of A and B;
355: /// ALPHA(1:K) = 1,
356: /// BETA(1:K) = 0,
357: /// and if M-K-L .GE. 0,
358: /// ALPHA(K+1:K+L) = diag(C),
359: /// BETA(K+1:K+L) = diag(S),
360: /// or if M-K-L .LT. 0,
361: /// ALPHA(K+1:M)= C, ALPHA(M+1:K+L)= 0
362: /// BETA(K+1:M) = S, BETA(M+1:K+L) = 1.
363: /// Furthermore, if K+L .LT. N,
364: /// ALPHA(K+L+1:N) = 0 and
365: /// BETA(K+L+1:N) = 0.
366: ///</param>
367: /// <param name="U">
368: /// (input/output) DOUBLE PRECISION array, dimension (LDU,M)
369: /// On entry, if JOBU = 'U', U must contain a matrix U1 (usually
370: /// the orthogonal matrix returned by DGGSVP).
371: /// On exit,
372: /// if JOBU = 'I', U contains the orthogonal matrix U;
373: /// if JOBU = 'U', U contains the product U1*U.
374: /// If JOBU = 'N', U is not referenced.
375: ///</param>
376: /// <param name="LDU">
377: /// (input) INTEGER
378: /// The leading dimension of the array U. LDU .GE. max(1,M) if
379: /// JOBU = 'U'; LDU .GE. 1 otherwise.
380: ///</param>
381: /// <param name="V">
382: /// (input/output) DOUBLE PRECISION array, dimension (LDV,P)
383: /// On entry, if JOBV = 'V', V must contain a matrix V1 (usually
384: /// the orthogonal matrix returned by DGGSVP).
385: /// On exit,
386: /// if JOBV = 'I', V contains the orthogonal matrix V;
387: /// if JOBV = 'V', V contains the product V1*V.
388: /// If JOBV = 'N', V is not referenced.
389: ///</param>
390: /// <param name="LDV">
391: /// (input) INTEGER
392: /// The leading dimension of the array V. LDV .GE. max(1,P) if
393: /// JOBV = 'V'; LDV .GE. 1 otherwise.
394: ///</param>
395: /// <param name="Q">
396: /// (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
397: /// On entry, if JOBQ = 'Q', Q must contain a matrix Q1 (usually
398: /// the orthogonal matrix returned by DGGSVP).
399: /// On exit,
400: /// if JOBQ = 'I', Q contains the orthogonal matrix Q;
401: /// if JOBQ = 'Q', Q contains the product Q1*Q.
402: /// If JOBQ = 'N', Q is not referenced.
403: ///</param>
404: /// <param name="LDQ">
405: /// (input) INTEGER
406: /// The leading dimension of the array Q. LDQ .GE. max(1,N) if
407: /// JOBQ = 'Q'; LDQ .GE. 1 otherwise.
408: ///</param>
409: /// <param name="WORK">
410: /// (workspace) DOUBLE PRECISION array, dimension (2*N)
411: ///</param>
412: /// <param name="NCYCLE">
413: /// (output) INTEGER
414: /// The number of cycles required for convergence.
415: ///</param>
416: /// <param name="INFO">
417: /// (output) INTEGER
418: /// = 0: successful exit
419: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
420: /// = 1: the procedure does not converge after MAXIT cycles.
421: ///</param>
422: public void Run(string JOBU, string JOBV, string JOBQ, int M, int P, int N
423: , int K, int L, ref double[] A, int offset_a, int LDA, ref double[] B, int offset_b, int LDB
424: , double TOLA, double TOLB, ref double[] ALPHA, int offset_alpha, ref double[] BETA, int offset_beta, ref double[] U, int offset_u, int LDU
425: , ref double[] V, int offset_v, int LDV, ref double[] Q, int offset_q, int LDQ, ref double[] WORK, int offset_work, ref int NCYCLE
426: , ref int INFO)
427: {
428:
429: #region Array Index Correction
430:
431: int o_a = -1 - LDA + offset_a; int o_b = -1 - LDB + offset_b; int o_alpha = -1 + offset_alpha;
432: int o_beta = -1 + offset_beta; int o_u = -1 - LDU + offset_u; int o_v = -1 - LDV + offset_v;
433: int o_q = -1 - LDQ + offset_q; int o_work = -1 + offset_work;
434:
435: #endregion
436:
437:
438: #region Strings
439:
440: JOBU = JOBU.Substring(0, 1); JOBV = JOBV.Substring(0, 1); JOBQ = JOBQ.Substring(0, 1);
441:
442: #endregion
443:
444:
445: #region Prolog
446:
447: // *
448: // * -- LAPACK routine (version 3.1) --
449: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
450: // * November 2006
451: // *
452: // * .. Scalar Arguments ..
453: // * ..
454: // * .. Array Arguments ..
455: // * ..
456: // *
457: // * Purpose
458: // * =======
459: // *
460: // * DTGSJA computes the generalized singular value decomposition (GSVD)
461: // * of two real upper triangular (or trapezoidal) matrices A and B.
462: // *
463: // * On entry, it is assumed that matrices A and B have the following
464: // * forms, which may be obtained by the preprocessing subroutine DGGSVP
465: // * from a general M-by-N matrix A and P-by-N matrix B:
466: // *
467: // * N-K-L K L
468: // * A = K ( 0 A12 A13 ) if M-K-L >= 0;
469: // * L ( 0 0 A23 )
470: // * M-K-L ( 0 0 0 )
471: // *
472: // * N-K-L K L
473: // * A = K ( 0 A12 A13 ) if M-K-L < 0;
474: // * M-K ( 0 0 A23 )
475: // *
476: // * N-K-L K L
477: // * B = L ( 0 0 B13 )
478: // * P-L ( 0 0 0 )
479: // *
480: // * where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular
481: // * upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0,
482: // * otherwise A23 is (M-K)-by-L upper trapezoidal.
483: // *
484: // * On exit,
485: // *
486: // * U'*A*Q = D1*( 0 R ), V'*B*Q = D2*( 0 R ),
487: // *
488: // * where U, V and Q are orthogonal matrices, Z' denotes the transpose
489: // * of Z, R is a nonsingular upper triangular matrix, and D1 and D2 are
490: // * ``diagonal'' matrices, which are of the following structures:
491: // *
492: // * If M-K-L >= 0,
493: // *
494: // * K L
495: // * D1 = K ( I 0 )
496: // * L ( 0 C )
497: // * M-K-L ( 0 0 )
498: // *
499: // * K L
500: // * D2 = L ( 0 S )
501: // * P-L ( 0 0 )
502: // *
503: // * N-K-L K L
504: // * ( 0 R ) = K ( 0 R11 R12 ) K
505: // * L ( 0 0 R22 ) L
506: // *
507: // * where
508: // *
509: // * C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
510: // * S = diag( BETA(K+1), ... , BETA(K+L) ),
511: // * C**2 + S**2 = I.
512: // *
513: // * R is stored in A(1:K+L,N-K-L+1:N) on exit.
514: // *
515: // * If M-K-L < 0,
516: // *
517: // * K M-K K+L-M
518: // * D1 = K ( I 0 0 )
519: // * M-K ( 0 C 0 )
520: // *
521: // * K M-K K+L-M
522: // * D2 = M-K ( 0 S 0 )
523: // * K+L-M ( 0 0 I )
524: // * P-L ( 0 0 0 )
525: // *
526: // * N-K-L K M-K K+L-M
527: // * ( 0 R ) = K ( 0 R11 R12 R13 )
528: // * M-K ( 0 0 R22 R23 )
529: // * K+L-M ( 0 0 0 R33 )
530: // *
531: // * where
532: // * C = diag( ALPHA(K+1), ... , ALPHA(M) ),
533: // * S = diag( BETA(K+1), ... , BETA(M) ),
534: // * C**2 + S**2 = I.
535: // *
536: // * R = ( R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N) and R33 is stored
537: // * ( 0 R22 R23 )
538: // * in B(M-K+1:L,N+M-K-L+1:N) on exit.
539: // *
540: // * The computation of the orthogonal transformation matrices U, V or Q
541: // * is optional. These matrices may either be formed explicitly, or they
542: // * may be postmultiplied into input matrices U1, V1, or Q1.
543: // *
544: // * Arguments
545: // * =========
546: // *
547: // * JOBU (input) CHARACTER*1
548: // * = 'U': U must contain an orthogonal matrix U1 on entry, and
549: // * the product U1*U is returned;
550: // * = 'I': U is initialized to the unit matrix, and the
551: // * orthogonal matrix U is returned;
552: // * = 'N': U is not computed.
553: // *
554: // * JOBV (input) CHARACTER*1
555: // * = 'V': V must contain an orthogonal matrix V1 on entry, and
556: // * the product V1*V is returned;
557: // * = 'I': V is initialized to the unit matrix, and the
558: // * orthogonal matrix V is returned;
559: // * = 'N': V is not computed.
560: // *
561: // * JOBQ (input) CHARACTER*1
562: // * = 'Q': Q must contain an orthogonal matrix Q1 on entry, and
563: // * the product Q1*Q is returned;
564: // * = 'I': Q is initialized to the unit matrix, and the
565: // * orthogonal matrix Q is returned;
566: // * = 'N': Q is not computed.
567: // *
568: // * M (input) INTEGER
569: // * The number of rows of the matrix A. M >= 0.
570: // *
571: // * P (input) INTEGER
572: // * The number of rows of the matrix B. P >= 0.
573: // *
574: // * N (input) INTEGER
575: // * The number of columns of the matrices A and B. N >= 0.
576: // *
577: // * K (input) INTEGER
578: // * L (input) INTEGER
579: // * K and L specify the subblocks in the input matrices A and B:
580: // * A23 = A(K+1:MIN(K+L,M),N-L+1:N) and B13 = B(1:L,N-L+1:N)
581: // * of A and B, whose GSVD is going to be computed by DTGSJA.
582: // * See Further details.
583: // *
584: // * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
585: // * On entry, the M-by-N matrix A.
586: // * On exit, A(N-K+1:N,1:MIN(K+L,M) ) contains the triangular
587: // * matrix R or part of R. See Purpose for details.
588: // *
589: // * LDA (input) INTEGER
590: // * The leading dimension of the array A. LDA >= max(1,M).
591: // *
592: // * B (input/output) DOUBLE PRECISION array, dimension (LDB,N)
593: // * On entry, the P-by-N matrix B.
594: // * On exit, if necessary, B(M-K+1:L,N+M-K-L+1:N) contains
595: // * a part of R. See Purpose for details.
596: // *
597: // * LDB (input) INTEGER
598: // * The leading dimension of the array B. LDB >= max(1,P).
599: // *
600: // * TOLA (input) DOUBLE PRECISION
601: // * TOLB (input) DOUBLE PRECISION
602: // * TOLA and TOLB are the convergence criteria for the Jacobi-
603: // * Kogbetliantz iteration procedure. Generally, they are the
604: // * same as used in the preprocessing step, say
605: // * TOLA = max(M,N)*norm(A)*MAZHEPS,
606: // * TOLB = max(P,N)*norm(B)*MAZHEPS.
607: // *
608: // * ALPHA (output) DOUBLE PRECISION array, dimension (N)
609: // * BETA (output) DOUBLE PRECISION array, dimension (N)
610: // * On exit, ALPHA and BETA contain the generalized singular
611: // * value pairs of A and B;
612: // * ALPHA(1:K) = 1,
613: // * BETA(1:K) = 0,
614: // * and if M-K-L >= 0,
615: // * ALPHA(K+1:K+L) = diag(C),
616: // * BETA(K+1:K+L) = diag(S),
617: // * or if M-K-L < 0,
618: // * ALPHA(K+1:M)= C, ALPHA(M+1:K+L)= 0
619: // * BETA(K+1:M) = S, BETA(M+1:K+L) = 1.
620: // * Furthermore, if K+L < N,
621: // * ALPHA(K+L+1:N) = 0 and
622: // * BETA(K+L+1:N) = 0.
623: // *
624: // * U (input/output) DOUBLE PRECISION array, dimension (LDU,M)
625: // * On entry, if JOBU = 'U', U must contain a matrix U1 (usually
626: // * the orthogonal matrix returned by DGGSVP).
627: // * On exit,
628: // * if JOBU = 'I', U contains the orthogonal matrix U;
629: // * if JOBU = 'U', U contains the product U1*U.
630: // * If JOBU = 'N', U is not referenced.
631: // *
632: // * LDU (input) INTEGER
633: // * The leading dimension of the array U. LDU >= max(1,M) if
634: // * JOBU = 'U'; LDU >= 1 otherwise.
635: // *
636: // * V (input/output) DOUBLE PRECISION array, dimension (LDV,P)
637: // * On entry, if JOBV = 'V', V must contain a matrix V1 (usually
638: // * the orthogonal matrix returned by DGGSVP).
639: // * On exit,
640: // * if JOBV = 'I', V contains the orthogonal matrix V;
641: // * if JOBV = 'V', V contains the product V1*V.
642: // * If JOBV = 'N', V is not referenced.
643: // *
644: // * LDV (input) INTEGER
645: // * The leading dimension of the array V. LDV >= max(1,P) if
646: // * JOBV = 'V'; LDV >= 1 otherwise.
647: // *
648: // * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
649: // * On entry, if JOBQ = 'Q', Q must contain a matrix Q1 (usually
650: // * the orthogonal matrix returned by DGGSVP).
651: // * On exit,
652: // * if JOBQ = 'I', Q contains the orthogonal matrix Q;
653: // * if JOBQ = 'Q', Q contains the product Q1*Q.
654: // * If JOBQ = 'N', Q is not referenced.
655: // *
656: // * LDQ (input) INTEGER
657: // * The leading dimension of the array Q. LDQ >= max(1,N) if
658: // * JOBQ = 'Q'; LDQ >= 1 otherwise.
659: // *
660: // * WORK (workspace) DOUBLE PRECISION array, dimension (2*N)
661: // *
662: // * NCYCLE (output) INTEGER
663: // * The number of cycles required for convergence.
664: // *
665: // * INFO (output) INTEGER
666: // * = 0: successful exit
667: // * < 0: if INFO = -i, the i-th argument had an illegal value.
668: // * = 1: the procedure does not converge after MAXIT cycles.
669: // *
670: // * Internal Parameters
671: // * ===================
672: // *
673: // * MAXIT INTEGER
674: // * MAXIT specifies the total loops that the iterative procedure
675: // * may take. If after MAXIT cycles, the routine fails to
676: // * converge, we return INFO = 1.
677: // *
678: // * Further Details
679: // * ===============
680: // *
681: // * DTGSJA essentially uses a variant of Kogbetliantz algorithm to reduce
682: // * min(L,M-K)-by-L triangular (or trapezoidal) matrix A23 and L-by-L
683: // * matrix B13 to the form:
684: // *
685: // * U1'*A13*Q1 = C1*R1; V1'*B13*Q1 = S1*R1,
686: // *
687: // * where U1, V1 and Q1 are orthogonal matrix, and Z' is the transpose
688: // * of Z. C1 and S1 are diagonal matrices satisfying
689: // *
690: // * C1**2 + S1**2 = I,
691: // *
692: // * and R1 is an L-by-L nonsingular upper triangular matrix.
693: // *
694: // * =====================================================================
695: // *
696: // * .. Parameters ..
697: // * ..
698: // * .. Local Scalars ..
699: // *
700: // * ..
701: // * .. External Functions ..
702: // * ..
703: // * .. External Subroutines ..
704: // * ..
705: // * .. Intrinsic Functions ..
706: // INTRINSIC ABS, MAX, MIN;
707: // * ..
708: // * .. Executable Statements ..
709: // *
710: // * Decode and test the input parameters
711: // *
712:
713: #endregion
714:
715:
716: #region Body
717:
718: INITU = this._lsame.Run(JOBU, "I");
719: WANTU = INITU || this._lsame.Run(JOBU, "U");
720: // *
721: INITV = this._lsame.Run(JOBV, "I");
722: WANTV = INITV || this._lsame.Run(JOBV, "V");
723: // *
724: INITQ = this._lsame.Run(JOBQ, "I");
725: WANTQ = INITQ || this._lsame.Run(JOBQ, "Q");
726: // *
727: INFO = 0;
728: if (!(INITU || WANTU || this._lsame.Run(JOBU, "N")))
729: {
730: INFO = - 1;
731: }
732: else
733: {
734: if (!(INITV || WANTV || this._lsame.Run(JOBV, "N")))
735: {
736: INFO = - 2;
737: }
738: else
739: {
740: if (!(INITQ || WANTQ || this._lsame.Run(JOBQ, "N")))
741: {
742: INFO = - 3;
743: }
744: else
745: {
746: if (M < 0)
747: {
748: INFO = - 4;
749: }
750: else
751: {
752: if (P < 0)
753: {
754: INFO = - 5;
755: }
756: else
757: {
758: if (N < 0)
759: {
760: INFO = - 6;
761: }
762: else
763: {
764: if (LDA < Math.Max(1, M))
765: {
766: INFO = - 10;
767: }
768: else
769: {
770: if (LDB < Math.Max(1, P))
771: {
772: INFO = - 12;
773: }
774: else
775: {
776: if (LDU < 1 || (WANTU && LDU < M))
777: {
778: INFO = - 18;
779: }
780: else
781: {
782: if (LDV < 1 || (WANTV && LDV < P))
783: {
784: INFO = - 20;
785: }
786: else
787: {
788: if (LDQ < 1 || (WANTQ && LDQ < N))
789: {
790: INFO = - 22;
791: }
792: }
793: }
794: }
795: }
796: }
797: }
798: }
799: }
800: }
801: }
802: if (INFO != 0)
803: {
804: this._xerbla.Run("DTGSJA", - INFO);
805: return;
806: }
807: // *
808: // * Initialize U, V and Q, if necessary
809: // *
810: if (INITU)
811: {
812: this._dlaset.Run("Full", M, M, ZERO, ONE, ref U, offset_u
813: , LDU);
814: }
815: if (INITV)
816: {
817: this._dlaset.Run("Full", P, P, ZERO, ONE, ref V, offset_v
818: , LDV);
819: }
820: if (INITQ)
821: {
822: this._dlaset.Run("Full", N, N, ZERO, ONE, ref Q, offset_q
823: , LDQ);
824: }
825: // *
826: // * Loop until convergence
827: // *
828: UPPER = false;
829: for (KCYCLE = 1; KCYCLE <= MAXIT; KCYCLE++)
830: {
831: // *
832: UPPER = !UPPER;
833: // *
834: for (I = 1; I <= L - 1; I++)
835: {
836: for (J = I + 1; J <= L; J++)
837: {
838: // *
839: A1 = ZERO;
840: A2 = ZERO;
841: A3 = ZERO;
842: if (K + I <= M) A1 = A[K + I+(N - L + I) * LDA + o_a];
843: if (K + J <= M) A3 = A[K + J+(N - L + J) * LDA + o_a];
844: // *
845: B1 = B[I+(N - L + I) * LDB + o_b];
846: B3 = B[J+(N - L + J) * LDB + o_b];
847: // *
848: if (UPPER)
849: {
850: if (K + I <= M) A2 = A[K + I+(N - L + J) * LDA + o_a];
851: B2 = B[I+(N - L + J) * LDB + o_b];
852: }
853: else
854: {
855: if (K + J <= M) A2 = A[K + J+(N - L + I) * LDA + o_a];
856: B2 = B[J+(N - L + I) * LDB + o_b];
857: }
858: // *
859: this._dlags2.Run(UPPER, A1, A2, A3, B1, B2
860: , B3, ref CSU, ref SNU, ref CSV, ref SNV, ref CSQ
861: , ref SNQ);
862: // *
863: // * Update (K+I)-th and (K+J)-th rows of matrix A: U'*A
864: // *
865: if (K + J <= M)
866: {
867: this._drot.Run(L, ref A, K + J+(N - L + 1) * LDA + o_a, LDA, ref A, K + I+(N - L + 1) * LDA + o_a, LDA, CSU
868: , SNU);
869: }
870: // *
871: // * Update I-th and J-th rows of matrix B: V'*B
872: // *
873: this._drot.Run(L, ref B, J+(N - L + 1) * LDB + o_b, LDB, ref B, I+(N - L + 1) * LDB + o_b, LDB, CSV
874: , SNV);
875: // *
876: // * Update (N-L+I)-th and (N-L+J)-th columns of matrices
877: // * A and B: A*Q and B*Q
878: // *
879: this._drot.Run(Math.Min(K + L, M), ref A, 1+(N - L + J) * LDA + o_a, 1, ref A, 1+(N - L + I) * LDA + o_a, 1, CSQ
880: , SNQ);
881: // *
882: this._drot.Run(L, ref B, 1+(N - L + J) * LDB + o_b, 1, ref B, 1+(N - L + I) * LDB + o_b, 1, CSQ
883: , SNQ);
884: // *
885: if (UPPER)
886: {
887: if (K + I <= M) A[K + I+(N - L + J) * LDA + o_a] = ZERO;
888: B[I+(N - L + J) * LDB + o_b] = ZERO;
889: }
890: else
891: {
892: if (K + J <= M) A[K + J+(N - L + I) * LDA + o_a] = ZERO;
893: B[J+(N - L + I) * LDB + o_b] = ZERO;
894: }
895: // *
896: // * Update orthogonal matrices U, V, Q, if desired.
897: // *
898: if (WANTU && K + J <= M)
899: {
900: this._drot.Run(M, ref U, 1+(K + J) * LDU + o_u, 1, ref U, 1+(K + I) * LDU + o_u, 1, CSU
901: , SNU);
902: }
903: // *
904: if (WANTV)
905: {
906: this._drot.Run(P, ref V, 1+J * LDV + o_v, 1, ref V, 1+I * LDV + o_v, 1, CSV
907: , SNV);
908: }
909: // *
910: if (WANTQ)
911: {
912: this._drot.Run(N, ref Q, 1+(N - L + J) * LDQ + o_q, 1, ref Q, 1+(N - L + I) * LDQ + o_q, 1, CSQ
913: , SNQ);
914: }
915: // *
916: }
917: }
918: // *
919: if (!UPPER)
920: {
921: // *
922: // * The matrices A13 and B13 were lower triangular at the start
923: // * of the cycle, and are now upper triangular.
924: // *
925: // * Convergence test: test the parallelism of the corresponding
926: // * rows of A and B.
927: // *
928: ERROR = ZERO;
929: for (I = 1; I <= Math.Min(L, M - K); I++)
930: {
931: this._dcopy.Run(L - I + 1, A, K + I+(N - L + I) * LDA + o_a, LDA, ref WORK, offset_work, 1);
932: this._dcopy.Run(L - I + 1, B, I+(N - L + I) * LDB + o_b, LDB, ref WORK, L + 1 + o_work, 1);
933: this._dlapll.Run(L - I + 1, ref WORK, offset_work, 1, ref WORK, L + 1 + o_work, 1, ref SSMIN);
934: ERROR = Math.Max(ERROR, SSMIN);
935: }
936: // *
937: if (Math.Abs(ERROR) <= Math.Min(TOLA, TOLB)) goto LABEL50;
938: }
939: // *
940: // * End of cycle loop
941: // *
942: }
943: // *
944: // * The algorithm has not converged after MAXIT cycles.
945: // *
946: INFO = 1;
947: goto LABEL100;
948: // *
949: LABEL50:;
950: // *
951: // * If ERROR <= MIN(TOLA,TOLB), then the algorithm has converged.
952: // * Compute the generalized singular value pairs (ALPHA, BETA), and
953: // * set the triangular matrix R to array A.
954: // *
955: for (I = 1; I <= K; I++)
956: {
957: ALPHA[I + o_alpha] = ONE;
958: BETA[I + o_beta] = ZERO;
959: }
960: // *
961: for (I = 1; I <= Math.Min(L, M - K); I++)
962: {
963: // *
964: A1 = A[K + I+(N - L + I) * LDA + o_a];
965: B1 = B[I+(N - L + I) * LDB + o_b];
966: // *
967: if (A1 != ZERO)
968: {
969: GAMMA = B1 / A1;
970: // *
971: // * change sign if necessary
972: // *
973: if (GAMMA < ZERO)
974: {
975: this._dscal.Run(L - I + 1, - ONE, ref B, I+(N - L + I) * LDB + o_b, LDB);
976: if (WANTV) this._dscal.Run(P, - ONE, ref V, 1+I * LDV + o_v, 1);
977: }
978: // *
979: this._dlartg.Run(Math.Abs(GAMMA), ONE, ref BETA[K + I + o_beta], ref ALPHA[K + I + o_alpha], ref RWK);
980: // *
981: if (ALPHA[K + I + o_alpha] >= BETA[K + I + o_beta])
982: {
983: this._dscal.Run(L - I + 1, ONE / ALPHA[K + I + o_alpha], ref A, K + I+(N - L + I) * LDA + o_a, LDA);
984: }
985: else
986: {
987: this._dscal.Run(L - I + 1, ONE / BETA[K + I + o_beta], ref B, I+(N - L + I) * LDB + o_b, LDB);
988: this._dcopy.Run(L - I + 1, B, I+(N - L + I) * LDB + o_b, LDB, ref A, K + I+(N - L + I) * LDA + o_a, LDA);
989: }
990: // *
991: }
992: else
993: {
994: // *
995: ALPHA[K + I + o_alpha] = ZERO;
996: BETA[K + I + o_beta] = ONE;
997: this._dcopy.Run(L - I + 1, B, I+(N - L + I) * LDB + o_b, LDB, ref A, K + I+(N - L + I) * LDA + o_a, LDA);
998: // *
999: }
1000: // *
1001: }
1002: // *
1003: // * Post-assignment
1004: // *
1005: for (I = M + 1; I <= K + L; I++)
1006: {
1007: ALPHA[I + o_alpha] = ZERO;
1008: BETA[I + o_beta] = ONE;
1009: }
1010: // *
1011: if (K + L < N)
1012: {
1013: for (I = K + L + 1; I <= N; I++)
1014: {
1015: ALPHA[I + o_alpha] = ZERO;
1016: BETA[I + o_beta] = ZERO;
1017: }
1018: }
1019: // *
1020: LABEL100:;
1021: NCYCLE = KCYCLE;
1022: return;
1023: // *
1024: // * End of DTGSJA
1025: // *
1026:
1027: #endregion
1028:
1029: }
1030: }
1031: }