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 driver routine (version 3.1) --
21: /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
22: /// November 2006
23: /// Purpose
24: /// =======
25: ///
26: /// DGGSVD computes the generalized singular value decomposition (GSVD)
27: /// of an M-by-N real matrix A and P-by-N real matrix B:
28: ///
29: /// U'*A*Q = D1*( 0 R ), V'*B*Q = D2*( 0 R )
30: ///
31: /// where U, V and Q are orthogonal matrices, and Z' is the transpose
32: /// of Z. Let K+L = the effective numerical rank of the matrix (A',B')',
33: /// then R is a K+L-by-K+L nonsingular upper triangular matrix, D1 and
34: /// D2 are M-by-(K+L) and P-by-(K+L) "diagonal" matrices and of the
35: /// following structures, respectively:
36: ///
37: /// If M-K-L .GE. 0,
38: ///
39: /// K L
40: /// D1 = K ( I 0 )
41: /// L ( 0 C )
42: /// M-K-L ( 0 0 )
43: ///
44: /// K L
45: /// D2 = L ( 0 S )
46: /// P-L ( 0 0 )
47: ///
48: /// N-K-L K L
49: /// ( 0 R ) = K ( 0 R11 R12 )
50: /// L ( 0 0 R22 )
51: ///
52: /// where
53: ///
54: /// C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
55: /// S = diag( BETA(K+1), ... , BETA(K+L) ),
56: /// C**2 + S**2 = I.
57: ///
58: /// R is stored in A(1:K+L,N-K-L+1:N) on exit.
59: ///
60: /// If M-K-L .LT. 0,
61: ///
62: /// K M-K K+L-M
63: /// D1 = K ( I 0 0 )
64: /// M-K ( 0 C 0 )
65: ///
66: /// K M-K K+L-M
67: /// D2 = M-K ( 0 S 0 )
68: /// K+L-M ( 0 0 I )
69: /// P-L ( 0 0 0 )
70: ///
71: /// N-K-L K M-K K+L-M
72: /// ( 0 R ) = K ( 0 R11 R12 R13 )
73: /// M-K ( 0 0 R22 R23 )
74: /// K+L-M ( 0 0 0 R33 )
75: ///
76: /// where
77: ///
78: /// C = diag( ALPHA(K+1), ... , ALPHA(M) ),
79: /// S = diag( BETA(K+1), ... , BETA(M) ),
80: /// C**2 + S**2 = I.
81: ///
82: /// (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored
83: /// ( 0 R22 R23 )
84: /// in B(M-K+1:L,N+M-K-L+1:N) on exit.
85: ///
86: /// The routine computes C, S, R, and optionally the orthogonal
87: /// transformation matrices U, V and Q.
88: ///
89: /// In particular, if B is an N-by-N nonsingular matrix, then the GSVD of
90: /// A and B implicitly gives the SVD of A*inv(B):
91: /// A*inv(B) = U*(D1*inv(D2))*V'.
92: /// If ( A',B')' has orthonormal columns, then the GSVD of A and B is
93: /// also equal to the CS decomposition of A and B. Furthermore, the GSVD
94: /// can be used to derive the solution of the eigenvalue problem:
95: /// A'*A x = lambda* B'*B x.
96: /// In some literature, the GSVD of A and B is presented in the form
97: /// U'*A*X = ( 0 D1 ), V'*B*X = ( 0 D2 )
98: /// where U and V are orthogonal and X is nonsingular, D1 and D2 are
99: /// ``diagonal''. The former GSVD form can be converted to the latter
100: /// form by taking the nonsingular matrix X as
101: ///
102: /// X = Q*( I 0 )
103: /// ( 0 inv(R) ).
104: ///
105: ///</summary>
106: public class DGGSVD
107: {
108:
109:
110: #region Dependencies
111:
112: LSAME _lsame; DLAMCH _dlamch; DLANGE _dlange; DCOPY _dcopy; DGGSVP _dggsvp; DTGSJA _dtgsja; XERBLA _xerbla;
113:
114: #endregion
115:
116:
117: #region Fields
118:
119: bool WANTQ = false; bool WANTU = false; bool WANTV = false; int I = 0; int IBND = 0; int ISUB = 0; int J = 0;
120: int NCYCLE = 0;double ANORM = 0; double BNORM = 0; double SMAX = 0; double TEMP = 0; double TOLA = 0; double TOLB = 0;
121: double ULP = 0;double UNFL = 0;
122:
123: #endregion
124:
125: public DGGSVD(LSAME lsame, DLAMCH dlamch, DLANGE dlange, DCOPY dcopy, DGGSVP dggsvp, DTGSJA dtgsja, XERBLA xerbla)
126: {
127:
128:
129: #region Set Dependencies
130:
131: this._lsame = lsame; this._dlamch = dlamch; this._dlange = dlange; this._dcopy = dcopy; this._dggsvp = dggsvp;
132: this._dtgsja = dtgsja;this._xerbla = xerbla;
133:
134: #endregion
135:
136: }
137:
138: public DGGSVD()
139: {
140:
141:
142: #region Dependencies (Initialization)
143:
144: LSAME lsame = new LSAME();
145: DLAMC3 dlamc3 = new DLAMC3();
146: DLASSQ dlassq = new DLASSQ();
147: DCOPY dcopy = new DCOPY();
148: XERBLA xerbla = new XERBLA();
149: DLAPY2 dlapy2 = new DLAPY2();
150: DNRM2 dnrm2 = new DNRM2();
151: DSCAL dscal = new DSCAL();
152: DSWAP dswap = new DSWAP();
153: IDAMAX idamax = new IDAMAX();
154: DLAPMT dlapmt = new DLAPMT();
155: DDOT ddot = new DDOT();
156: DAXPY daxpy = new DAXPY();
157: DLAS2 dlas2 = new DLAS2();
158: DROT drot = new DROT();
159: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
160: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
161: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
162: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
163: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
164: DLANGE dlange = new DLANGE(dlassq, lsame);
165: DGEMV dgemv = new DGEMV(lsame, xerbla);
166: DGER dger = new DGER(xerbla);
167: DLARF dlarf = new DLARF(dgemv, dger, lsame);
168: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
169: DGEQR2 dgeqr2 = new DGEQR2(dlarf, dlarfg, xerbla);
170: DORM2R dorm2r = new DORM2R(lsame, dlarf, xerbla);
171: DGEQPF dgeqpf = new DGEQPF(dgeqr2, dlarf, dlarfg, dorm2r, dswap, xerbla, idamax, dlamch, dnrm2);
172: DGERQ2 dgerq2 = new DGERQ2(dlarf, dlarfg, xerbla);
173: DLACPY dlacpy = new DLACPY(lsame);
174: DLASET dlaset = new DLASET(lsame);
175: DORG2R dorg2r = new DORG2R(dlarf, dscal, xerbla);
176: DORMR2 dormr2 = new DORMR2(lsame, dlarf, xerbla);
177: DGGSVP dggsvp = new DGGSVP(lsame, dgeqpf, dgeqr2, dgerq2, dlacpy, dlapmt, dlaset, dorg2r, dorm2r, dormr2
178: , xerbla);
179: DLARTG dlartg = new DLARTG(dlamch);
180: DLASV2 dlasv2 = new DLASV2(dlamch);
181: DLAGS2 dlags2 = new DLAGS2(dlartg, dlasv2);
182: DLAPLL dlapll = new DLAPLL(ddot, daxpy, dlarfg, dlas2);
183: DTGSJA dtgsja = new DTGSJA(lsame, dcopy, dlags2, dlapll, dlartg, dlaset, drot, dscal, xerbla);
184:
185: #endregion
186:
187:
188: #region Set Dependencies
189:
190: this._lsame = lsame; this._dlamch = dlamch; this._dlange = dlange; this._dcopy = dcopy; this._dggsvp = dggsvp;
191: this._dtgsja = dtgsja;this._xerbla = xerbla;
192:
193: #endregion
194:
195: }
196: /// <summary>
197: /// Purpose
198: /// =======
199: ///
200: /// DGGSVD computes the generalized singular value decomposition (GSVD)
201: /// of an M-by-N real matrix A and P-by-N real matrix B:
202: ///
203: /// U'*A*Q = D1*( 0 R ), V'*B*Q = D2*( 0 R )
204: ///
205: /// where U, V and Q are orthogonal matrices, and Z' is the transpose
206: /// of Z. Let K+L = the effective numerical rank of the matrix (A',B')',
207: /// then R is a K+L-by-K+L nonsingular upper triangular matrix, D1 and
208: /// D2 are M-by-(K+L) and P-by-(K+L) "diagonal" matrices and of the
209: /// following structures, respectively:
210: ///
211: /// If M-K-L .GE. 0,
212: ///
213: /// K L
214: /// D1 = K ( I 0 )
215: /// L ( 0 C )
216: /// M-K-L ( 0 0 )
217: ///
218: /// K L
219: /// D2 = L ( 0 S )
220: /// P-L ( 0 0 )
221: ///
222: /// N-K-L K L
223: /// ( 0 R ) = K ( 0 R11 R12 )
224: /// L ( 0 0 R22 )
225: ///
226: /// where
227: ///
228: /// C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
229: /// S = diag( BETA(K+1), ... , BETA(K+L) ),
230: /// C**2 + S**2 = I.
231: ///
232: /// R is stored in A(1:K+L,N-K-L+1:N) on exit.
233: ///
234: /// If M-K-L .LT. 0,
235: ///
236: /// K M-K K+L-M
237: /// D1 = K ( I 0 0 )
238: /// M-K ( 0 C 0 )
239: ///
240: /// K M-K K+L-M
241: /// D2 = M-K ( 0 S 0 )
242: /// K+L-M ( 0 0 I )
243: /// P-L ( 0 0 0 )
244: ///
245: /// N-K-L K M-K K+L-M
246: /// ( 0 R ) = K ( 0 R11 R12 R13 )
247: /// M-K ( 0 0 R22 R23 )
248: /// K+L-M ( 0 0 0 R33 )
249: ///
250: /// where
251: ///
252: /// C = diag( ALPHA(K+1), ... , ALPHA(M) ),
253: /// S = diag( BETA(K+1), ... , BETA(M) ),
254: /// C**2 + S**2 = I.
255: ///
256: /// (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored
257: /// ( 0 R22 R23 )
258: /// in B(M-K+1:L,N+M-K-L+1:N) on exit.
259: ///
260: /// The routine computes C, S, R, and optionally the orthogonal
261: /// transformation matrices U, V and Q.
262: ///
263: /// In particular, if B is an N-by-N nonsingular matrix, then the GSVD of
264: /// A and B implicitly gives the SVD of A*inv(B):
265: /// A*inv(B) = U*(D1*inv(D2))*V'.
266: /// If ( A',B')' has orthonormal columns, then the GSVD of A and B is
267: /// also equal to the CS decomposition of A and B. Furthermore, the GSVD
268: /// can be used to derive the solution of the eigenvalue problem:
269: /// A'*A x = lambda* B'*B x.
270: /// In some literature, the GSVD of A and B is presented in the form
271: /// U'*A*X = ( 0 D1 ), V'*B*X = ( 0 D2 )
272: /// where U and V are orthogonal and X is nonsingular, D1 and D2 are
273: /// ``diagonal''. The former GSVD form can be converted to the latter
274: /// form by taking the nonsingular matrix X as
275: ///
276: /// X = Q*( I 0 )
277: /// ( 0 inv(R) ).
278: ///
279: ///</summary>
280: /// <param name="JOBU">
281: /// (input) CHARACTER*1
282: /// = 'U': Orthogonal matrix U is computed;
283: /// = 'N': U is not computed.
284: ///</param>
285: /// <param name="JOBV">
286: /// (input) CHARACTER*1
287: /// = 'V': Orthogonal matrix V is computed;
288: /// = 'N': V is not computed.
289: ///</param>
290: /// <param name="JOBQ">
291: /// (input) CHARACTER*1
292: /// = 'Q': Orthogonal matrix Q is computed;
293: /// = 'N': Q is not computed.
294: ///</param>
295: /// <param name="M">
296: /// (input) INTEGER
297: /// The number of rows of the matrix A. M .GE. 0.
298: ///</param>
299: /// <param name="N">
300: /// (input) INTEGER
301: /// The number of columns of the matrices A and B. N .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="K">
308: /// L
309: ///</param>
310: /// <param name="L">
311: /// ( 0 C )
312: ///</param>
313: /// <param name="A">
314: /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
315: /// On entry, the M-by-N matrix A.
316: /// On exit, A contains the triangular matrix R, or part of R.
317: /// See Purpose for details.
318: ///</param>
319: /// <param name="LDA">
320: /// (input) INTEGER
321: /// The leading dimension of the array A. LDA .GE. max(1,M).
322: ///</param>
323: /// <param name="B">
324: /// (input/output) DOUBLE PRECISION array, dimension (LDB,N)
325: /// On entry, the P-by-N matrix B.
326: /// On exit, B contains the triangular matrix R if M-K-L .LT. 0.
327: /// See Purpose for details.
328: ///</param>
329: /// <param name="LDB">
330: /// (input) INTEGER
331: /// The leading dimension of the array B. LDB .GE. max(1,P).
332: ///</param>
333: /// <param name="ALPHA">
334: /// (output) DOUBLE PRECISION array, dimension (N)
335: ///</param>
336: /// <param name="BETA">
337: /// (output) DOUBLE PRECISION array, dimension (N)
338: /// On exit, ALPHA and BETA contain the generalized singular
339: /// value pairs of A and B;
340: /// ALPHA(1:K) = 1,
341: /// BETA(1:K) = 0,
342: /// and if M-K-L .GE. 0,
343: /// ALPHA(K+1:K+L) = C,
344: /// BETA(K+1:K+L) = S,
345: /// or if M-K-L .LT. 0,
346: /// ALPHA(K+1:M)=C, ALPHA(M+1:K+L)=0
347: /// BETA(K+1:M) =S, BETA(M+1:K+L) =1
348: /// and
349: /// ALPHA(K+L+1:N) = 0
350: /// BETA(K+L+1:N) = 0
351: ///</param>
352: /// <param name="U">
353: /// (output) DOUBLE PRECISION array, dimension (LDU,M)
354: /// If JOBU = 'U', U contains the M-by-M orthogonal matrix U.
355: /// If JOBU = 'N', U is not referenced.
356: ///</param>
357: /// <param name="LDU">
358: /// (input) INTEGER
359: /// The leading dimension of the array U. LDU .GE. max(1,M) if
360: /// JOBU = 'U'; LDU .GE. 1 otherwise.
361: ///</param>
362: /// <param name="V">
363: /// (output) DOUBLE PRECISION array, dimension (LDV,P)
364: /// If JOBV = 'V', V contains the P-by-P orthogonal matrix V.
365: /// If JOBV = 'N', V is not referenced.
366: ///</param>
367: /// <param name="LDV">
368: /// (input) INTEGER
369: /// The leading dimension of the array V. LDV .GE. max(1,P) if
370: /// JOBV = 'V'; LDV .GE. 1 otherwise.
371: ///</param>
372: /// <param name="Q">
373: /// (output) DOUBLE PRECISION array, dimension (LDQ,N)
374: /// If JOBQ = 'Q', Q contains the N-by-N orthogonal matrix Q.
375: /// If JOBQ = 'N', Q is not referenced.
376: ///</param>
377: /// <param name="LDQ">
378: /// (input) INTEGER
379: /// The leading dimension of the array Q. LDQ .GE. max(1,N) if
380: /// JOBQ = 'Q'; LDQ .GE. 1 otherwise.
381: ///</param>
382: /// <param name="WORK">
383: /// (workspace) DOUBLE PRECISION array,
384: /// dimension (max(3*N,M,P)+N)
385: ///</param>
386: /// <param name="IWORK">
387: /// (workspace/output) INTEGER array, dimension (N)
388: /// On exit, IWORK stores the sorting information. More
389: /// precisely, the following loop will sort ALPHA
390: /// for I = K+1, min(M,K+L)
391: /// swap ALPHA(I) and ALPHA(IWORK(I))
392: /// endfor
393: /// such that ALPHA(1) .GE. ALPHA(2) .GE. ... .GE. ALPHA(N).
394: ///</param>
395: /// <param name="INFO">
396: /// (output) INTEGER
397: /// = 0: successful exit
398: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
399: /// .GT. 0: if INFO = 1, the Jacobi-type procedure failed to
400: /// converge. For further details, see subroutine DTGSJA.
401: ///</param>
402: public void Run(string JOBU, string JOBV, string JOBQ, int M, int N, int P
403: , ref int K, ref int L, ref double[] A, int offset_a, int LDA, ref double[] B, int offset_b, int LDB
404: , ref double[] ALPHA, int offset_alpha, ref double[] BETA, int offset_beta, ref double[] U, int offset_u, int LDU, ref double[] V, int offset_v, int LDV
405: , ref double[] Q, int offset_q, int LDQ, ref double[] WORK, int offset_work, ref int[] IWORK, int offset_iwork, ref int INFO)
406: {
407:
408: #region Array Index Correction
409:
410: int o_a = -1 - LDA + offset_a; int o_b = -1 - LDB + offset_b; int o_alpha = -1 + offset_alpha;
411: int o_beta = -1 + offset_beta; int o_u = -1 - LDU + offset_u; int o_v = -1 - LDV + offset_v;
412: int o_q = -1 - LDQ + offset_q; int o_work = -1 + offset_work; int o_iwork = -1 + offset_iwork;
413:
414: #endregion
415:
416:
417: #region Strings
418:
419: JOBU = JOBU.Substring(0, 1); JOBV = JOBV.Substring(0, 1); JOBQ = JOBQ.Substring(0, 1);
420:
421: #endregion
422:
423:
424: #region Prolog
425:
426: // *
427: // * -- LAPACK driver routine (version 3.1) --
428: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
429: // * November 2006
430: // *
431: // * .. Scalar Arguments ..
432: // * ..
433: // * .. Array Arguments ..
434: // * ..
435: // *
436: // * Purpose
437: // * =======
438: // *
439: // * DGGSVD computes the generalized singular value decomposition (GSVD)
440: // * of an M-by-N real matrix A and P-by-N real matrix B:
441: // *
442: // * U'*A*Q = D1*( 0 R ), V'*B*Q = D2*( 0 R )
443: // *
444: // * where U, V and Q are orthogonal matrices, and Z' is the transpose
445: // * of Z. Let K+L = the effective numerical rank of the matrix (A',B')',
446: // * then R is a K+L-by-K+L nonsingular upper triangular matrix, D1 and
447: // * D2 are M-by-(K+L) and P-by-(K+L) "diagonal" matrices and of the
448: // * following structures, respectively:
449: // *
450: // * If M-K-L >= 0,
451: // *
452: // * K L
453: // * D1 = K ( I 0 )
454: // * L ( 0 C )
455: // * M-K-L ( 0 0 )
456: // *
457: // * K L
458: // * D2 = L ( 0 S )
459: // * P-L ( 0 0 )
460: // *
461: // * N-K-L K L
462: // * ( 0 R ) = K ( 0 R11 R12 )
463: // * L ( 0 0 R22 )
464: // *
465: // * where
466: // *
467: // * C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
468: // * S = diag( BETA(K+1), ... , BETA(K+L) ),
469: // * C**2 + S**2 = I.
470: // *
471: // * R is stored in A(1:K+L,N-K-L+1:N) on exit.
472: // *
473: // * If M-K-L < 0,
474: // *
475: // * K M-K K+L-M
476: // * D1 = K ( I 0 0 )
477: // * M-K ( 0 C 0 )
478: // *
479: // * K M-K K+L-M
480: // * D2 = M-K ( 0 S 0 )
481: // * K+L-M ( 0 0 I )
482: // * P-L ( 0 0 0 )
483: // *
484: // * N-K-L K M-K K+L-M
485: // * ( 0 R ) = K ( 0 R11 R12 R13 )
486: // * M-K ( 0 0 R22 R23 )
487: // * K+L-M ( 0 0 0 R33 )
488: // *
489: // * where
490: // *
491: // * C = diag( ALPHA(K+1), ... , ALPHA(M) ),
492: // * S = diag( BETA(K+1), ... , BETA(M) ),
493: // * C**2 + S**2 = I.
494: // *
495: // * (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored
496: // * ( 0 R22 R23 )
497: // * in B(M-K+1:L,N+M-K-L+1:N) on exit.
498: // *
499: // * The routine computes C, S, R, and optionally the orthogonal
500: // * transformation matrices U, V and Q.
501: // *
502: // * In particular, if B is an N-by-N nonsingular matrix, then the GSVD of
503: // * A and B implicitly gives the SVD of A*inv(B):
504: // * A*inv(B) = U*(D1*inv(D2))*V'.
505: // * If ( A',B')' has orthonormal columns, then the GSVD of A and B is
506: // * also equal to the CS decomposition of A and B. Furthermore, the GSVD
507: // * can be used to derive the solution of the eigenvalue problem:
508: // * A'*A x = lambda* B'*B x.
509: // * In some literature, the GSVD of A and B is presented in the form
510: // * U'*A*X = ( 0 D1 ), V'*B*X = ( 0 D2 )
511: // * where U and V are orthogonal and X is nonsingular, D1 and D2 are
512: // * ``diagonal''. The former GSVD form can be converted to the latter
513: // * form by taking the nonsingular matrix X as
514: // *
515: // * X = Q*( I 0 )
516: // * ( 0 inv(R) ).
517: // *
518: // * Arguments
519: // * =========
520: // *
521: // * JOBU (input) CHARACTER*1
522: // * = 'U': Orthogonal matrix U is computed;
523: // * = 'N': U is not computed.
524: // *
525: // * JOBV (input) CHARACTER*1
526: // * = 'V': Orthogonal matrix V is computed;
527: // * = 'N': V is not computed.
528: // *
529: // * JOBQ (input) CHARACTER*1
530: // * = 'Q': Orthogonal matrix Q is computed;
531: // * = 'N': Q is not computed.
532: // *
533: // * M (input) INTEGER
534: // * The number of rows of the matrix A. M >= 0.
535: // *
536: // * N (input) INTEGER
537: // * The number of columns of the matrices A and B. N >= 0.
538: // *
539: // * P (input) INTEGER
540: // * The number of rows of the matrix B. P >= 0.
541: // *
542: // * K (output) INTEGER
543: // * L (output) INTEGER
544: // * On exit, K and L specify the dimension of the subblocks
545: // * described in the Purpose section.
546: // * K + L = effective numerical rank of (A',B')'.
547: // *
548: // * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
549: // * On entry, the M-by-N matrix A.
550: // * On exit, A contains the triangular matrix R, or part of R.
551: // * See Purpose for details.
552: // *
553: // * LDA (input) INTEGER
554: // * The leading dimension of the array A. LDA >= max(1,M).
555: // *
556: // * B (input/output) DOUBLE PRECISION array, dimension (LDB,N)
557: // * On entry, the P-by-N matrix B.
558: // * On exit, B contains the triangular matrix R if M-K-L < 0.
559: // * See Purpose for details.
560: // *
561: // * LDB (input) INTEGER
562: // * The leading dimension of the array B. LDB >= max(1,P).
563: // *
564: // * ALPHA (output) DOUBLE PRECISION array, dimension (N)
565: // * BETA (output) DOUBLE PRECISION array, dimension (N)
566: // * On exit, ALPHA and BETA contain the generalized singular
567: // * value pairs of A and B;
568: // * ALPHA(1:K) = 1,
569: // * BETA(1:K) = 0,
570: // * and if M-K-L >= 0,
571: // * ALPHA(K+1:K+L) = C,
572: // * BETA(K+1:K+L) = S,
573: // * or if M-K-L < 0,
574: // * ALPHA(K+1:M)=C, ALPHA(M+1:K+L)=0
575: // * BETA(K+1:M) =S, BETA(M+1:K+L) =1
576: // * and
577: // * ALPHA(K+L+1:N) = 0
578: // * BETA(K+L+1:N) = 0
579: // *
580: // * U (output) DOUBLE PRECISION array, dimension (LDU,M)
581: // * If JOBU = 'U', U contains the M-by-M orthogonal matrix U.
582: // * If JOBU = 'N', U is not referenced.
583: // *
584: // * LDU (input) INTEGER
585: // * The leading dimension of the array U. LDU >= max(1,M) if
586: // * JOBU = 'U'; LDU >= 1 otherwise.
587: // *
588: // * V (output) DOUBLE PRECISION array, dimension (LDV,P)
589: // * If JOBV = 'V', V contains the P-by-P orthogonal matrix V.
590: // * If JOBV = 'N', V is not referenced.
591: // *
592: // * LDV (input) INTEGER
593: // * The leading dimension of the array V. LDV >= max(1,P) if
594: // * JOBV = 'V'; LDV >= 1 otherwise.
595: // *
596: // * Q (output) DOUBLE PRECISION array, dimension (LDQ,N)
597: // * If JOBQ = 'Q', Q contains the N-by-N orthogonal matrix Q.
598: // * If JOBQ = 'N', Q is not referenced.
599: // *
600: // * LDQ (input) INTEGER
601: // * The leading dimension of the array Q. LDQ >= max(1,N) if
602: // * JOBQ = 'Q'; LDQ >= 1 otherwise.
603: // *
604: // * WORK (workspace) DOUBLE PRECISION array,
605: // * dimension (max(3*N,M,P)+N)
606: // *
607: // * IWORK (workspace/output) INTEGER array, dimension (N)
608: // * On exit, IWORK stores the sorting information. More
609: // * precisely, the following loop will sort ALPHA
610: // * for I = K+1, min(M,K+L)
611: // * swap ALPHA(I) and ALPHA(IWORK(I))
612: // * endfor
613: // * such that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(N).
614: // *
615: // * INFO (output) INTEGER
616: // * = 0: successful exit
617: // * < 0: if INFO = -i, the i-th argument had an illegal value.
618: // * > 0: if INFO = 1, the Jacobi-type procedure failed to
619: // * converge. For further details, see subroutine DTGSJA.
620: // *
621: // * Internal Parameters
622: // * ===================
623: // *
624: // * TOLA DOUBLE PRECISION
625: // * TOLB DOUBLE PRECISION
626: // * TOLA and TOLB are the thresholds to determine the effective
627: // * rank of (A',B')'. Generally, they are set to
628: // * TOLA = MAX(M,N)*norm(A)*MAZHEPS,
629: // * TOLB = MAX(P,N)*norm(B)*MAZHEPS.
630: // * The size of TOLA and TOLB may affect the size of backward
631: // * errors of the decomposition.
632: // *
633: // * Further Details
634: // * ===============
635: // *
636: // * 2-96 Based on modifications by
637: // * Ming Gu and Huan Ren, Computer Science Division, University of
638: // * California at Berkeley, USA
639: // *
640: // * =====================================================================
641: // *
642: // * .. Local Scalars ..
643: // * ..
644: // * .. External Functions ..
645: // * ..
646: // * .. External Subroutines ..
647: // * ..
648: // * .. Intrinsic Functions ..
649: // INTRINSIC MAX, MIN;
650: // * ..
651: // * .. Executable Statements ..
652: // *
653: // * Test the input parameters
654: // *
655:
656: #endregion
657:
658:
659: #region Body
660:
661: WANTU = this._lsame.Run(JOBU, "U");
662: WANTV = this._lsame.Run(JOBV, "V");
663: WANTQ = this._lsame.Run(JOBQ, "Q");
664: // *
665: INFO = 0;
666: if (!(WANTU || this._lsame.Run(JOBU, "N")))
667: {
668: INFO = - 1;
669: }
670: else
671: {
672: if (!(WANTV || this._lsame.Run(JOBV, "N")))
673: {
674: INFO = - 2;
675: }
676: else
677: {
678: if (!(WANTQ || this._lsame.Run(JOBQ, "N")))
679: {
680: INFO = - 3;
681: }
682: else
683: {
684: if (M < 0)
685: {
686: INFO = - 4;
687: }
688: else
689: {
690: if (N < 0)
691: {
692: INFO = - 5;
693: }
694: else
695: {
696: if (P < 0)
697: {
698: INFO = - 6;
699: }
700: else
701: {
702: if (LDA < Math.Max(1, M))
703: {
704: INFO = - 10;
705: }
706: else
707: {
708: if (LDB < Math.Max(1, P))
709: {
710: INFO = - 12;
711: }
712: else
713: {
714: if (LDU < 1 || (WANTU && LDU < M))
715: {
716: INFO = - 16;
717: }
718: else
719: {
720: if (LDV < 1 || (WANTV && LDV < P))
721: {
722: INFO = - 18;
723: }
724: else
725: {
726: if (LDQ < 1 || (WANTQ && LDQ < N))
727: {
728: INFO = - 20;
729: }
730: }
731: }
732: }
733: }
734: }
735: }
736: }
737: }
738: }
739: }
740: if (INFO != 0)
741: {
742: this._xerbla.Run("DGGSVD", - INFO);
743: return;
744: }
745: // *
746: // * Compute the Frobenius norm of matrices A and B
747: // *
748: ANORM = this._dlange.Run("1", M, N, A, offset_a, LDA, ref WORK, offset_work);
749: BNORM = this._dlange.Run("1", P, N, B, offset_b, LDB, ref WORK, offset_work);
750: // *
751: // * Get machine precision and set up threshold for determining
752: // * the effective numerical rank of the matrices A and B.
753: // *
754: ULP = this._dlamch.Run("Precision");
755: UNFL = this._dlamch.Run("Safe Minimum");
756: TOLA = Math.Max(M, N) * Math.Max(ANORM, UNFL) * ULP;
757: TOLB = Math.Max(P, N) * Math.Max(BNORM, UNFL) * ULP;
758: // *
759: // * Preprocessing
760: // *
761: this._dggsvp.Run(JOBU, JOBV, JOBQ, M, P, N
762: , ref A, offset_a, LDA, ref B, offset_b, LDB, TOLA, TOLB
763: , ref K, ref L, ref U, offset_u, LDU, ref V, offset_v, LDV
764: , ref Q, offset_q, LDQ, ref IWORK, offset_iwork, ref WORK, offset_work, ref WORK, N + 1 + o_work, ref INFO);
765: // *
766: // * Compute the GSVD of two upper "triangular" matrices
767: // *
768: this._dtgsja.Run(JOBU, JOBV, JOBQ, M, P, N
769: , K, L, ref A, offset_a, LDA, ref B, offset_b, LDB
770: , TOLA, TOLB, ref ALPHA, offset_alpha, ref BETA, offset_beta, ref U, offset_u, LDU
771: , ref V, offset_v, LDV, ref Q, offset_q, LDQ, ref WORK, offset_work, ref NCYCLE
772: , ref INFO);
773: // *
774: // * Sort the singular values and store the pivot indices in IWORK
775: // * Copy ALPHA to WORK, then sort ALPHA in WORK
776: // *
777: this._dcopy.Run(N, ALPHA, offset_alpha, 1, ref WORK, offset_work, 1);
778: IBND = Math.Min(L, M - K);
779: for (I = 1; I <= IBND; I++)
780: {
781: // *
782: // * Scan for largest ALPHA(K+I)
783: // *
784: ISUB = I;
785: SMAX = WORK[K + I + o_work];
786: for (J = I + 1; J <= IBND; J++)
787: {
788: TEMP = WORK[K + J + o_work];
789: if (TEMP > SMAX)
790: {
791: ISUB = J;
792: SMAX = TEMP;
793: }
794: }
795: if (ISUB != I)
796: {
797: WORK[K + ISUB + o_work] = WORK[K + I + o_work];
798: WORK[K + I + o_work] = SMAX;
799: IWORK[K + I + o_iwork] = K + ISUB;
800: }
801: else
802: {
803: IWORK[K + I + o_iwork] = K + I;
804: }
805: }
806: // *
807: return;
808: // *
809: // * End of DGGSVD
810: // *
811:
812: #endregion
813:
814: }
815: }
816: }