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: /// DGESDD computes the singular value decomposition (SVD) of a real
27: /// M-by-N matrix A, optionally computing the left and right singular
28: /// vectors. If singular vectors are desired, it uses a
29: /// divide-and-conquer algorithm.
30: ///
31: /// The SVD is written
32: ///
33: /// A = U * SIGMA * transpose(V)
34: ///
35: /// where SIGMA is an M-by-N matrix which is zero except for its
36: /// min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
37: /// V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA
38: /// are the singular values of A; they are real and non-negative, and
39: /// are returned in descending order. The first min(m,n) columns of
40: /// U and V are the left and right singular vectors of A.
41: ///
42: /// Note that the routine returns VT = V**T, not V.
43: ///
44: /// The divide and conquer algorithm makes very mild assumptions about
45: /// floating point arithmetic. It will work on machines with a guard
46: /// digit in add/subtract, or on those binary machines without guard
47: /// digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
48: /// Cray-2. It could conceivably fail on hexadecimal or decimal machines
49: /// without guard digits, but we know of none.
50: ///
51: ///</summary>
52: public class DGESDD
53: {
54:
55:
56: #region Dependencies
57:
58: DBDSDC _dbdsdc; DGEBRD _dgebrd; DGELQF _dgelqf; DGEMM _dgemm; DGEQRF _dgeqrf; DLACPY _dlacpy; DLASCL _dlascl;
59: DLASET _dlaset;DORGBR _dorgbr; DORGLQ _dorglq; DORGQR _dorgqr; DORMBR _dormbr; XERBLA _xerbla; DLAMCH _dlamch;
60: DLANGE _dlange;ILAENV _ilaenv; LSAME _lsame;
61:
62: #endregion
63:
64:
65: #region Fields
66:
67: const double ZERO = 0.0E0; const double ONE = 1.0E0; bool LQUERY = false; bool WNTQA = false; bool WNTQAS = false;
68: bool WNTQN = false;bool WNTQO = false; bool WNTQS = false; int BDSPAC = 0; int BLK = 0; int CHUNK = 0; int I = 0;
69: int IE = 0;int IERR = 0; int IL = 0; int IR = 0; int ISCL = 0; int ITAU = 0; int ITAUP = 0; int ITAUQ = 0; int IU = 0;
70: int IVT = 0;int LDWKVT = 0; int LDWRKL = 0; int LDWRKR = 0; int LDWRKU = 0; int MAXWRK = 0; int MINMN = 0; int MINWRK = 0;
71: int MNTHR = 0;int NWORK = 0; int WRKBL = 0; double ANRM = 0; double BIGNUM = 0; double EPS = 0; double SMLNUM = 0;
72: int[] IDUM = new int[1]; int offset_idum = 0; double[] DUM = new double[1]; int offset_dum = 0;
73:
74: #endregion
75:
76: public DGESDD(DBDSDC dbdsdc, DGEBRD dgebrd, DGELQF dgelqf, DGEMM dgemm, DGEQRF dgeqrf, DLACPY dlacpy, DLASCL dlascl, DLASET dlaset, DORGBR dorgbr, DORGLQ dorglq
77: , DORGQR dorgqr, DORMBR dormbr, XERBLA xerbla, DLAMCH dlamch, DLANGE dlange, ILAENV ilaenv, LSAME lsame)
78: {
79:
80:
81: #region Set Dependencies
82:
83: this._dbdsdc = dbdsdc; this._dgebrd = dgebrd; this._dgelqf = dgelqf; this._dgemm = dgemm; this._dgeqrf = dgeqrf;
84: this._dlacpy = dlacpy;this._dlascl = dlascl; this._dlaset = dlaset; this._dorgbr = dorgbr; this._dorglq = dorglq;
85: this._dorgqr = dorgqr;this._dormbr = dormbr; this._xerbla = xerbla; this._dlamch = dlamch; this._dlange = dlange;
86: this._ilaenv = ilaenv;this._lsame = lsame;
87:
88: #endregion
89:
90: }
91:
92: public DGESDD()
93: {
94:
95:
96: #region Dependencies (Initialization)
97:
98: LSAME lsame = new LSAME();
99: IEEECK ieeeck = new IEEECK();
100: IPARMQ iparmq = new IPARMQ();
101: DLAMC3 dlamc3 = new DLAMC3();
102: DLASSQ dlassq = new DLASSQ();
103: DCOPY dcopy = new DCOPY();
104: XERBLA xerbla = new XERBLA();
105: DLAMRG dlamrg = new DLAMRG();
106: DLAPY2 dlapy2 = new DLAPY2();
107: DROT drot = new DROT();
108: DNRM2 dnrm2 = new DNRM2();
109: DLASD5 dlasd5 = new DLASD5();
110: DLAS2 dlas2 = new DLAS2();
111: DLASQ5 dlasq5 = new DLASQ5();
112: DLAZQ4 dlazq4 = new DLAZQ4();
113: DSCAL dscal = new DSCAL();
114: DSWAP dswap = new DSWAP();
115: DLASDT dlasdt = new DLASDT();
116: DDOT ddot = new DDOT();
117: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
118: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
119: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
120: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
121: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
122: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
123: DLANST dlanst = new DLANST(lsame, dlassq);
124: DLARTG dlartg = new DLARTG(dlamch);
125: DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
126: DLACPY dlacpy = new DLACPY(lsame);
127: DLASET dlaset = new DLASET(lsame);
128: DLASD2 dlasd2 = new DLASD2(dlamch, dlapy2, dcopy, dlacpy, dlamrg, dlaset, drot, xerbla);
129: DGEMM dgemm = new DGEMM(lsame, xerbla);
130: DLAED6 dlaed6 = new DLAED6(dlamch);
131: DLASD4 dlasd4 = new DLASD4(dlaed6, dlasd5, dlamch);
132: DLASD3 dlasd3 = new DLASD3(dlamc3, dnrm2, dcopy, dgemm, dlacpy, dlascl, dlasd4, xerbla);
133: DLASD1 dlasd1 = new DLASD1(dlamrg, dlascl, dlasd2, dlasd3, xerbla);
134: DLASQ6 dlasq6 = new DLASQ6(dlamch);
135: DLAZQ3 dlazq3 = new DLAZQ3(dlasq5, dlasq6, dlazq4, dlamch);
136: DLASRT dlasrt = new DLASRT(lsame, xerbla);
137: DLASQ2 dlasq2 = new DLASQ2(dlazq3, dlasrt, xerbla, dlamch, ilaenv);
138: DLASQ1 dlasq1 = new DLASQ1(dcopy, dlas2, dlascl, dlasq2, dlasrt, xerbla, dlamch);
139: DLASR dlasr = new DLASR(lsame, xerbla);
140: DLASV2 dlasv2 = new DLASV2(dlamch);
141: DBDSQR dbdsqr = new DBDSQR(lsame, dlamch, dlartg, dlas2, dlasq1, dlasr, dlasv2, drot, dscal, dswap
142: , xerbla);
143: DLASDQ dlasdq = new DLASDQ(dbdsqr, dlartg, dlasr, dswap, xerbla, lsame);
144: DLASD0 dlasd0 = new DLASD0(dlasd1, dlasdq, dlasdt, xerbla);
145: DLASD7 dlasd7 = new DLASD7(dcopy, dlamrg, drot, xerbla, dlamch, dlapy2);
146: DLASD8 dlasd8 = new DLASD8(dcopy, dlascl, dlasd4, dlaset, xerbla, ddot, dlamc3, dnrm2);
147: DLASD6 dlasd6 = new DLASD6(dcopy, dlamrg, dlascl, dlasd7, dlasd8, xerbla);
148: DLASDA dlasda = new DLASDA(dcopy, dlasd6, dlasdq, dlasdt, dlaset, xerbla);
149: DBDSDC dbdsdc = new DBDSDC(lsame, ilaenv, dlamch, dlanst, dcopy, dlartg, dlascl, dlasd0, dlasda, dlasdq
150: , dlaset, dlasr, dswap, xerbla);
151: DGEMV dgemv = new DGEMV(lsame, xerbla);
152: DGER dger = new DGER(xerbla);
153: DLARF dlarf = new DLARF(dgemv, dger, lsame);
154: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
155: DGEBD2 dgebd2 = new DGEBD2(dlarf, dlarfg, xerbla);
156: DLABRD dlabrd = new DLABRD(dgemv, dlarfg, dscal);
157: DGEBRD dgebrd = new DGEBRD(dgebd2, dgemm, dlabrd, xerbla, ilaenv);
158: DGELQ2 dgelq2 = new DGELQ2(dlarf, dlarfg, xerbla);
159: DTRMM dtrmm = new DTRMM(lsame, xerbla);
160: DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
161: DTRMV dtrmv = new DTRMV(lsame, xerbla);
162: DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);
163: DGELQF dgelqf = new DGELQF(dgelq2, dlarfb, dlarft, xerbla, ilaenv);
164: DGEQR2 dgeqr2 = new DGEQR2(dlarf, dlarfg, xerbla);
165: DGEQRF dgeqrf = new DGEQRF(dgeqr2, dlarfb, dlarft, xerbla, ilaenv);
166: DORGL2 dorgl2 = new DORGL2(dlarf, dscal, xerbla);
167: DORGLQ dorglq = new DORGLQ(dlarfb, dlarft, dorgl2, xerbla, ilaenv);
168: DORG2R dorg2r = new DORG2R(dlarf, dscal, xerbla);
169: DORGQR dorgqr = new DORGQR(dlarfb, dlarft, dorg2r, xerbla, ilaenv);
170: DORGBR dorgbr = new DORGBR(lsame, ilaenv, dorglq, dorgqr, xerbla);
171: DORML2 dorml2 = new DORML2(lsame, dlarf, xerbla);
172: DORMLQ dormlq = new DORMLQ(lsame, ilaenv, dlarfb, dlarft, dorml2, xerbla);
173: DORM2R dorm2r = new DORM2R(lsame, dlarf, xerbla);
174: DORMQR dormqr = new DORMQR(lsame, ilaenv, dlarfb, dlarft, dorm2r, xerbla);
175: DORMBR dormbr = new DORMBR(lsame, ilaenv, dormlq, dormqr, xerbla);
176: DLANGE dlange = new DLANGE(dlassq, lsame);
177:
178: #endregion
179:
180:
181: #region Set Dependencies
182:
183: this._dbdsdc = dbdsdc; this._dgebrd = dgebrd; this._dgelqf = dgelqf; this._dgemm = dgemm; this._dgeqrf = dgeqrf;
184: this._dlacpy = dlacpy;this._dlascl = dlascl; this._dlaset = dlaset; this._dorgbr = dorgbr; this._dorglq = dorglq;
185: this._dorgqr = dorgqr;this._dormbr = dormbr; this._xerbla = xerbla; this._dlamch = dlamch; this._dlange = dlange;
186: this._ilaenv = ilaenv;this._lsame = lsame;
187:
188: #endregion
189:
190: }
191: /// <summary>
192: /// Purpose
193: /// =======
194: ///
195: /// DGESDD computes the singular value decomposition (SVD) of a real
196: /// M-by-N matrix A, optionally computing the left and right singular
197: /// vectors. If singular vectors are desired, it uses a
198: /// divide-and-conquer algorithm.
199: ///
200: /// The SVD is written
201: ///
202: /// A = U * SIGMA * transpose(V)
203: ///
204: /// where SIGMA is an M-by-N matrix which is zero except for its
205: /// min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
206: /// V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA
207: /// are the singular values of A; they are real and non-negative, and
208: /// are returned in descending order. The first min(m,n) columns of
209: /// U and V are the left and right singular vectors of A.
210: ///
211: /// Note that the routine returns VT = V**T, not V.
212: ///
213: /// The divide and conquer algorithm makes very mild assumptions about
214: /// floating point arithmetic. It will work on machines with a guard
215: /// digit in add/subtract, or on those binary machines without guard
216: /// digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
217: /// Cray-2. It could conceivably fail on hexadecimal or decimal machines
218: /// without guard digits, but we know of none.
219: ///
220: ///</summary>
221: /// <param name="JOBZ">
222: /// (input) CHARACTER*1
223: /// Specifies options for computing all or part of the matrix U:
224: /// = 'A': all M columns of U and all N rows of V**T are
225: /// returned in the arrays U and VT;
226: /// = 'S': the first min(M,N) columns of U and the first
227: /// min(M,N) rows of V**T are returned in the arrays U
228: /// and VT;
229: /// = 'O': If M .GE. N, the first N columns of U are overwritten
230: /// on the array A and all rows of V**T are returned in
231: /// the array VT;
232: /// otherwise, all columns of U are returned in the
233: /// array U and the first M rows of V**T are overwritten
234: /// in the array A;
235: /// = 'N': no columns of U or rows of V**T are computed.
236: ///</param>
237: /// <param name="M">
238: /// (input) INTEGER
239: /// The number of rows of the input matrix A. M .GE. 0.
240: ///</param>
241: /// <param name="N">
242: /// (input) INTEGER
243: /// The number of columns of the input matrix A. N .GE. 0.
244: ///</param>
245: /// <param name="A">
246: /// = U * SIGMA * transpose(V)
247: ///</param>
248: /// <param name="LDA">
249: /// (input) INTEGER
250: /// The leading dimension of the array A. LDA .GE. max(1,M).
251: ///</param>
252: /// <param name="S">
253: /// (output) DOUBLE PRECISION array, dimension (min(M,N))
254: /// The singular values of A, sorted so that S(i) .GE. S(i+1).
255: ///</param>
256: /// <param name="U">
257: /// (output) DOUBLE PRECISION array, dimension (LDU,UCOL)
258: /// UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M .LT. N;
259: /// UCOL = min(M,N) if JOBZ = 'S'.
260: /// If JOBZ = 'A' or JOBZ = 'O' and M .LT. N, U contains the M-by-M
261: /// orthogonal matrix U;
262: /// if JOBZ = 'S', U contains the first min(M,N) columns of U
263: /// (the left singular vectors, stored columnwise);
264: /// if JOBZ = 'O' and M .GE. N, or JOBZ = 'N', U is not referenced.
265: ///</param>
266: /// <param name="LDU">
267: /// (input) INTEGER
268: /// The leading dimension of the array U. LDU .GE. 1; if
269: /// JOBZ = 'S' or 'A' or JOBZ = 'O' and M .LT. N, LDU .GE. M.
270: ///</param>
271: /// <param name="VT">
272: /// (output) DOUBLE PRECISION array, dimension (LDVT,N)
273: /// If JOBZ = 'A' or JOBZ = 'O' and M .GE. N, VT contains the
274: /// N-by-N orthogonal matrix V**T;
275: /// if JOBZ = 'S', VT contains the first min(M,N) rows of
276: /// V**T (the right singular vectors, stored rowwise);
277: /// if JOBZ = 'O' and M .LT. N, or JOBZ = 'N', VT is not referenced.
278: ///</param>
279: /// <param name="LDVT">
280: /// (input) INTEGER
281: /// The leading dimension of the array VT. LDVT .GE. 1; if
282: /// JOBZ = 'A' or JOBZ = 'O' and M .GE. N, LDVT .GE. N;
283: /// if JOBZ = 'S', LDVT .GE. min(M,N).
284: ///</param>
285: /// <param name="WORK">
286: /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
287: /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
288: ///</param>
289: /// <param name="LWORK">
290: /// (input) INTEGER
291: /// The dimension of the array WORK. LWORK .GE. 1.
292: /// If JOBZ = 'N',
293: /// LWORK .GE. 3*min(M,N) + max(max(M,N),7*min(M,N)).
294: /// If JOBZ = 'O',
295: /// LWORK .GE. 3*min(M,N)*min(M,N) +
296: /// max(max(M,N),5*min(M,N)*min(M,N)+4*min(M,N)).
297: /// If JOBZ = 'S' or 'A'
298: /// LWORK .GE. 3*min(M,N)*min(M,N) +
299: /// max(max(M,N),4*min(M,N)*min(M,N)+4*min(M,N)).
300: /// For good performance, LWORK should generally be larger.
301: /// If LWORK = -1 but other input arguments are legal, WORK(1)
302: /// returns the optimal LWORK.
303: ///</param>
304: /// <param name="IWORK">
305: /// (workspace) INTEGER array, dimension (8*min(M,N))
306: ///</param>
307: /// <param name="INFO">
308: /// (output) INTEGER
309: /// = 0: successful exit.
310: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
311: /// .GT. 0: DBDSDC did not converge, updating process failed.
312: ///</param>
313: public void Run(string JOBZ, int M, int N, ref double[] A, int offset_a, int LDA, ref double[] S, int offset_s
314: , ref double[] U, int offset_u, int LDU, ref double[] VT, int offset_vt, int LDVT, ref double[] WORK, int offset_work, int LWORK
315: , ref int[] IWORK, int offset_iwork, ref int INFO)
316: {
317:
318: #region Array Index Correction
319:
320: int o_a = -1 - LDA + offset_a; int o_s = -1 + offset_s; int o_u = -1 - LDU + offset_u;
321: int o_vt = -1 - LDVT + offset_vt; int o_work = -1 + offset_work; int o_iwork = -1 + offset_iwork;
322:
323: #endregion
324:
325:
326: #region Strings
327:
328: JOBZ = JOBZ.Substring(0, 1);
329:
330: #endregion
331:
332:
333: #region Prolog
334:
335: // *
336: // * -- LAPACK driver routine (version 3.1) --
337: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
338: // * November 2006
339: // *
340: // * .. Scalar Arguments ..
341: // * ..
342: // * .. Array Arguments ..
343: // * ..
344: // *
345: // * Purpose
346: // * =======
347: // *
348: // * DGESDD computes the singular value decomposition (SVD) of a real
349: // * M-by-N matrix A, optionally computing the left and right singular
350: // * vectors. If singular vectors are desired, it uses a
351: // * divide-and-conquer algorithm.
352: // *
353: // * The SVD is written
354: // *
355: // * A = U * SIGMA * transpose(V)
356: // *
357: // * where SIGMA is an M-by-N matrix which is zero except for its
358: // * min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
359: // * V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA
360: // * are the singular values of A; they are real and non-negative, and
361: // * are returned in descending order. The first min(m,n) columns of
362: // * U and V are the left and right singular vectors of A.
363: // *
364: // * Note that the routine returns VT = V**T, not V.
365: // *
366: // * The divide and conquer algorithm makes very mild assumptions about
367: // * floating point arithmetic. It will work on machines with a guard
368: // * digit in add/subtract, or on those binary machines without guard
369: // * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
370: // * Cray-2. It could conceivably fail on hexadecimal or decimal machines
371: // * without guard digits, but we know of none.
372: // *
373: // * Arguments
374: // * =========
375: // *
376: // * JOBZ (input) CHARACTER*1
377: // * Specifies options for computing all or part of the matrix U:
378: // * = 'A': all M columns of U and all N rows of V**T are
379: // * returned in the arrays U and VT;
380: // * = 'S': the first min(M,N) columns of U and the first
381: // * min(M,N) rows of V**T are returned in the arrays U
382: // * and VT;
383: // * = 'O': If M >= N, the first N columns of U are overwritten
384: // * on the array A and all rows of V**T are returned in
385: // * the array VT;
386: // * otherwise, all columns of U are returned in the
387: // * array U and the first M rows of V**T are overwritten
388: // * in the array A;
389: // * = 'N': no columns of U or rows of V**T are computed.
390: // *
391: // * M (input) INTEGER
392: // * The number of rows of the input matrix A. M >= 0.
393: // *
394: // * N (input) INTEGER
395: // * The number of columns of the input matrix A. N >= 0.
396: // *
397: // * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
398: // * On entry, the M-by-N matrix A.
399: // * On exit,
400: // * if JOBZ = 'O', A is overwritten with the first N columns
401: // * of U (the left singular vectors, stored
402: // * columnwise) if M >= N;
403: // * A is overwritten with the first M rows
404: // * of V**T (the right singular vectors, stored
405: // * rowwise) otherwise.
406: // * if JOBZ .ne. 'O', the contents of A are destroyed.
407: // *
408: // * LDA (input) INTEGER
409: // * The leading dimension of the array A. LDA >= max(1,M).
410: // *
411: // * S (output) DOUBLE PRECISION array, dimension (min(M,N))
412: // * The singular values of A, sorted so that S(i) >= S(i+1).
413: // *
414: // * U (output) DOUBLE PRECISION array, dimension (LDU,UCOL)
415: // * UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N;
416: // * UCOL = min(M,N) if JOBZ = 'S'.
417: // * If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M
418: // * orthogonal matrix U;
419: // * if JOBZ = 'S', U contains the first min(M,N) columns of U
420: // * (the left singular vectors, stored columnwise);
421: // * if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.
422: // *
423: // * LDU (input) INTEGER
424: // * The leading dimension of the array U. LDU >= 1; if
425: // * JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.
426: // *
427: // * VT (output) DOUBLE PRECISION array, dimension (LDVT,N)
428: // * If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the
429: // * N-by-N orthogonal matrix V**T;
430: // * if JOBZ = 'S', VT contains the first min(M,N) rows of
431: // * V**T (the right singular vectors, stored rowwise);
432: // * if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.
433: // *
434: // * LDVT (input) INTEGER
435: // * The leading dimension of the array VT. LDVT >= 1; if
436: // * JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;
437: // * if JOBZ = 'S', LDVT >= min(M,N).
438: // *
439: // * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
440: // * On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
441: // *
442: // * LWORK (input) INTEGER
443: // * The dimension of the array WORK. LWORK >= 1.
444: // * If JOBZ = 'N',
445: // * LWORK >= 3*min(M,N) + max(max(M,N),7*min(M,N)).
446: // * If JOBZ = 'O',
447: // * LWORK >= 3*min(M,N)*min(M,N) +
448: // * max(max(M,N),5*min(M,N)*min(M,N)+4*min(M,N)).
449: // * If JOBZ = 'S' or 'A'
450: // * LWORK >= 3*min(M,N)*min(M,N) +
451: // * max(max(M,N),4*min(M,N)*min(M,N)+4*min(M,N)).
452: // * For good performance, LWORK should generally be larger.
453: // * If LWORK = -1 but other input arguments are legal, WORK(1)
454: // * returns the optimal LWORK.
455: // *
456: // * IWORK (workspace) INTEGER array, dimension (8*min(M,N))
457: // *
458: // * INFO (output) INTEGER
459: // * = 0: successful exit.
460: // * < 0: if INFO = -i, the i-th argument had an illegal value.
461: // * > 0: DBDSDC did not converge, updating process failed.
462: // *
463: // * Further Details
464: // * ===============
465: // *
466: // * Based on contributions by
467: // * Ming Gu and Huan Ren, Computer Science Division, University of
468: // * California at Berkeley, USA
469: // *
470: // * =====================================================================
471: // *
472: // * .. Parameters ..
473: // * ..
474: // * .. Local Scalars ..
475: // * ..
476: // * .. Local Arrays ..
477: // * ..
478: // * .. External Subroutines ..
479: // * ..
480: // * .. External Functions ..
481: // * ..
482: // * .. Intrinsic Functions ..
483: // INTRINSIC INT, MAX, MIN, SQRT;
484: // * ..
485: // * .. Executable Statements ..
486: // *
487: // * Test the input arguments
488: // *
489:
490: #endregion
491:
492:
493: #region Body
494:
495: INFO = 0;
496: MINMN = Math.Min(M, N);
497: WNTQA = this._lsame.Run(JOBZ, "A");
498: WNTQS = this._lsame.Run(JOBZ, "S");
499: WNTQAS = WNTQA || WNTQS;
500: WNTQO = this._lsame.Run(JOBZ, "O");
501: WNTQN = this._lsame.Run(JOBZ, "N");
502: LQUERY = (LWORK == - 1);
503: // *
504: if (!(WNTQA || WNTQS || WNTQO || WNTQN))
505: {
506: INFO = - 1;
507: }
508: else
509: {
510: if (M < 0)
511: {
512: INFO = - 2;
513: }
514: else
515: {
516: if (N < 0)
517: {
518: INFO = - 3;
519: }
520: else
521: {
522: if (LDA < Math.Max(1, M))
523: {
524: INFO = - 5;
525: }
526: else
527: {
528: if (LDU < 1 || (WNTQAS && LDU < M) || (WNTQO && M < N && LDU < M))
529: {
530: INFO = - 8;
531: }
532: else
533: {
534: if (LDVT < 1 || (WNTQA && LDVT < N) || (WNTQS && LDVT < MINMN) || (WNTQO && M >= N && LDVT < N))
535: {
536: INFO = - 10;
537: }
538: }
539: }
540: }
541: }
542: }
543: // *
544: // * Compute workspace
545: // * (Note: Comments in the code beginning "Workspace:" describe the
546: // * minimal amount of workspace needed at that point in the code,
547: // * as well as the preferred amount for good performance.
548: // * NB refers to the optimal block size for the immediately
549: // * following subroutine, as returned by ILAENV.)
550: // *
551: if (INFO == 0)
552: {
553: MINWRK = 1;
554: MAXWRK = 1;
555: if (M >= N && MINMN > 0)
556: {
557: // *
558: // * Compute space needed for DBDSDC
559: // *
560: MNTHR = Convert.ToInt32(Math.Truncate(MINMN * 11.0E0 / 6.0E0));
561: if (WNTQN)
562: {
563: BDSPAC = 7 * N;
564: }
565: else
566: {
567: BDSPAC = 3 * N * N + 4 * N;
568: }
569: if (M >= MNTHR)
570: {
571: if (WNTQN)
572: {
573: // *
574: // * Path 1 (M much larger than N, JOBZ='N')
575: // *
576: WRKBL = N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N, - 1, - 1);
577: WRKBL = Math.Max(WRKBL, 3 * N + 2 * N * this._ilaenv.Run(1, "DGEBRD", " ", N, N, - 1, - 1));
578: MAXWRK = Math.Max(WRKBL, BDSPAC + N);
579: MINWRK = BDSPAC + N;
580: }
581: else
582: {
583: if (WNTQO)
584: {
585: // *
586: // * Path 2 (M much larger than N, JOBZ='O')
587: // *
588: WRKBL = N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N, - 1, - 1);
589: WRKBL = Math.Max(WRKBL, N + N * this._ilaenv.Run(1, "DORGQR", " ", M, N, N, - 1));
590: WRKBL = Math.Max(WRKBL, 3 * N + 2 * N * this._ilaenv.Run(1, "DGEBRD", " ", N, N, - 1, - 1));
591: WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORMBR", "QLN", N, N, N, - 1));
592: WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORMBR", "PRT", N, N, N, - 1));
593: WRKBL = Math.Max(WRKBL, BDSPAC + 3 * N);
594: MAXWRK = WRKBL + 2 * N * N;
595: MINWRK = BDSPAC + 2 * N * N + 3 * N;
596: }
597: else
598: {
599: if (WNTQS)
600: {
601: // *
602: // * Path 3 (M much larger than N, JOBZ='S')
603: // *
604: WRKBL = N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N, - 1, - 1);
605: WRKBL = Math.Max(WRKBL, N + N * this._ilaenv.Run(1, "DORGQR", " ", M, N, N, - 1));
606: WRKBL = Math.Max(WRKBL, 3 * N + 2 * N * this._ilaenv.Run(1, "DGEBRD", " ", N, N, - 1, - 1));
607: WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORMBR", "QLN", N, N, N, - 1));
608: WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORMBR", "PRT", N, N, N, - 1));
609: WRKBL = Math.Max(WRKBL, BDSPAC + 3 * N);
610: MAXWRK = WRKBL + N * N;
611: MINWRK = BDSPAC + N * N + 3 * N;
612: }
613: else
614: {
615: if (WNTQA)
616: {
617: // *
618: // * Path 4 (M much larger than N, JOBZ='A')
619: // *
620: WRKBL = N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N, - 1, - 1);
621: WRKBL = Math.Max(WRKBL, N + M * this._ilaenv.Run(1, "DORGQR", " ", M, M, N, - 1));
622: WRKBL = Math.Max(WRKBL, 3 * N + 2 * N * this._ilaenv.Run(1, "DGEBRD", " ", N, N, - 1, - 1));
623: WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORMBR", "QLN", N, N, N, - 1));
624: WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORMBR", "PRT", N, N, N, - 1));
625: WRKBL = Math.Max(WRKBL, BDSPAC + 3 * N);
626: MAXWRK = WRKBL + N * N;
627: MINWRK = BDSPAC + N * N + 3 * N;
628: }
629: }
630: }
631: }
632: }
633: else
634: {
635: // *
636: // * Path 5 (M at least N, but not much larger)
637: // *
638: WRKBL = 3 * N + (M + N) * this._ilaenv.Run(1, "DGEBRD", " ", M, N, - 1, - 1);
639: if (WNTQN)
640: {
641: MAXWRK = Math.Max(WRKBL, BDSPAC + 3 * N);
642: MINWRK = 3 * N + Math.Max(M, BDSPAC);
643: }
644: else
645: {
646: if (WNTQO)
647: {
648: WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORMBR", "QLN", M, N, N, - 1));
649: WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORMBR", "PRT", N, N, N, - 1));
650: WRKBL = Math.Max(WRKBL, BDSPAC + 3 * N);
651: MAXWRK = WRKBL + M * N;
652: MINWRK = 3 * N + Math.Max(M, N * N + BDSPAC);
653: }
654: else
655: {
656: if (WNTQS)
657: {
658: WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORMBR", "QLN", M, N, N, - 1));
659: WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORMBR", "PRT", N, N, N, - 1));
660: MAXWRK = Math.Max(WRKBL, BDSPAC + 3 * N);
661: MINWRK = 3 * N + Math.Max(M, BDSPAC);
662: }
663: else
664: {
665: if (WNTQA)
666: {
667: WRKBL = Math.Max(WRKBL, 3 * N + M * this._ilaenv.Run(1, "DORMBR", "QLN", M, M, N, - 1));
668: WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORMBR", "PRT", N, N, N, - 1));
669: MAXWRK = Math.Max(MAXWRK, BDSPAC + 3 * N);
670: MINWRK = 3 * N + Math.Max(M, BDSPAC);
671: }
672: }
673: }
674: }
675: }
676: }
677: else
678: {
679: if (MINMN > 0)
680: {
681: // *
682: // * Compute space needed for DBDSDC
683: // *
684: MNTHR = Convert.ToInt32(Math.Truncate(MINMN * 11.0E0 / 6.0E0));
685: if (WNTQN)
686: {
687: BDSPAC = 7 * M;
688: }
689: else
690: {
691: BDSPAC = 3 * M * M + 4 * M;
692: }
693: if (N >= MNTHR)
694: {
695: if (WNTQN)
696: {
697: // *
698: // * Path 1t (N much larger than M, JOBZ='N')
699: // *
700: WRKBL = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N, - 1, - 1);
701: WRKBL = Math.Max(WRKBL, 3 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M, - 1, - 1));
702: MAXWRK = Math.Max(WRKBL, BDSPAC + M);
703: MINWRK = BDSPAC + M;
704: }
705: else
706: {
707: if (WNTQO)
708: {
709: // *
710: // * Path 2t (N much larger than M, JOBZ='O')
711: // *
712: WRKBL = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N, - 1, - 1);
713: WRKBL = Math.Max(WRKBL, M + M * this._ilaenv.Run(1, "DORGLQ", " ", M, N, M, - 1));
714: WRKBL = Math.Max(WRKBL, 3 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M, - 1, - 1));
715: WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORMBR", "QLN", M, M, M, - 1));
716: WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORMBR", "PRT", M, M, M, - 1));
717: WRKBL = Math.Max(WRKBL, BDSPAC + 3 * M);
718: MAXWRK = WRKBL + 2 * M * M;
719: MINWRK = BDSPAC + 2 * M * M + 3 * M;
720: }
721: else
722: {
723: if (WNTQS)
724: {
725: // *
726: // * Path 3t (N much larger than M, JOBZ='S')
727: // *
728: WRKBL = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N, - 1, - 1);
729: WRKBL = Math.Max(WRKBL, M + M * this._ilaenv.Run(1, "DORGLQ", " ", M, N, M, - 1));
730: WRKBL = Math.Max(WRKBL, 3 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M, - 1, - 1));
731: WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORMBR", "QLN", M, M, M, - 1));
732: WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORMBR", "PRT", M, M, M, - 1));
733: WRKBL = Math.Max(WRKBL, BDSPAC + 3 * M);
734: MAXWRK = WRKBL + M * M;
735: MINWRK = BDSPAC + M * M + 3 * M;
736: }
737: else
738: {
739: if (WNTQA)
740: {
741: // *
742: // * Path 4t (N much larger than M, JOBZ='A')
743: // *
744: WRKBL = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N, - 1, - 1);
745: WRKBL = Math.Max(WRKBL, M + N * this._ilaenv.Run(1, "DORGLQ", " ", N, N, M, - 1));
746: WRKBL = Math.Max(WRKBL, 3 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M, - 1, - 1));
747: WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORMBR", "QLN", M, M, M, - 1));
748: WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORMBR", "PRT", M, M, M, - 1));
749: WRKBL = Math.Max(WRKBL, BDSPAC + 3 * M);
750: MAXWRK = WRKBL + M * M;
751: MINWRK = BDSPAC + M * M + 3 * M;
752: }
753: }
754: }
755: }
756: }
757: else
758: {
759: // *
760: // * Path 5t (N greater than M, but not much larger)
761: // *
762: WRKBL = 3 * M + (M + N) * this._ilaenv.Run(1, "DGEBRD", " ", M, N, - 1, - 1);
763: if (WNTQN)
764: {
765: MAXWRK = Math.Max(WRKBL, BDSPAC + 3 * M);
766: MINWRK = 3 * M + Math.Max(N, BDSPAC);
767: }
768: else
769: {
770: if (WNTQO)
771: {
772: WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORMBR", "QLN", M, M, N, - 1));
773: WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORMBR", "PRT", M, N, M, - 1));
774: WRKBL = Math.Max(WRKBL, BDSPAC + 3 * M);
775: MAXWRK = WRKBL + M * N;
776: MINWRK = 3 * M + Math.Max(N, M * M + BDSPAC);
777: }
778: else
779: {
780: if (WNTQS)
781: {
782: WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORMBR", "QLN", M, M, N, - 1));
783: WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORMBR", "PRT", M, N, M, - 1));
784: MAXWRK = Math.Max(WRKBL, BDSPAC + 3 * M);
785: MINWRK = 3 * M + Math.Max(N, BDSPAC);
786: }
787: else
788: {
789: if (WNTQA)
790: {
791: WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORMBR", "QLN", M, M, N, - 1));
792: WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORMBR", "PRT", N, N, M, - 1));
793: MAXWRK = Math.Max(WRKBL, BDSPAC + 3 * M);
794: MINWRK = 3 * M + Math.Max(N, BDSPAC);
795: }
796: }
797: }
798: }
799: }
800: }
801: }
802: MAXWRK = Math.Max(MAXWRK, MINWRK);
803: WORK[1 + o_work] = MAXWRK;
804: // *
805: if (LWORK < MINWRK && !LQUERY)
806: {
807: INFO = - 12;
808: }
809: }
810: // *
811: if (INFO != 0)
812: {
813: this._xerbla.Run("DGESDD", - INFO);
814: return;
815: }
816: else
817: {
818: if (LQUERY)
819: {
820: return;
821: }
822: }
823: // *
824: // * Quick return if possible
825: // *
826: if (M == 0 || N == 0)
827: {
828: return;
829: }
830: // *
831: // * Get machine constants
832: // *
833: EPS = this._dlamch.Run("P");
834: SMLNUM = Math.Sqrt(this._dlamch.Run("S")) / EPS;
835: BIGNUM = ONE / SMLNUM;
836: // *
837: // * Scale A if max element outside range [SMLNUM,BIGNUM]
838: // *
839: ANRM = this._dlange.Run("M", M, N, A, offset_a, LDA, ref DUM, offset_dum);
840: ISCL = 0;
841: if (ANRM > ZERO && ANRM < SMLNUM)
842: {
843: ISCL = 1;
844: this._dlascl.Run("G", 0, 0, ANRM, SMLNUM, M
845: , N, ref A, offset_a, LDA, ref IERR);
846: }
847: else
848: {
849: if (ANRM > BIGNUM)
850: {
851: ISCL = 1;
852: this._dlascl.Run("G", 0, 0, ANRM, BIGNUM, M
853: , N, ref A, offset_a, LDA, ref IERR);
854: }
855: }
856: // *
857: if (M >= N)
858: {
859: // *
860: // * A has at least as many rows as columns. If A has sufficiently
861: // * more rows than columns, first reduce using the QR
862: // * decomposition (if sufficient workspace available)
863: // *
864: if (M >= MNTHR)
865: {
866: // *
867: if (WNTQN)
868: {
869: // *
870: // * Path 1 (M much larger than N, JOBZ='N')
871: // * No singular vectors to be computed
872: // *
873: ITAU = 1;
874: NWORK = ITAU + N;
875: // *
876: // * Compute A=Q*R
877: // * (Workspace: need 2*N, prefer N+N*NB)
878: // *
879: this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, NWORK + o_work
880: , LWORK - NWORK + 1, ref IERR);
881: // *
882: // * Zero out below R
883: // *
884: this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref A, 2+1 * LDA + o_a
885: , LDA);
886: IE = 1;
887: ITAUQ = IE + N;
888: ITAUP = ITAUQ + N;
889: NWORK = ITAUP + N;
890: // *
891: // * Bidiagonalize R in A
892: // * (Workspace: need 4*N, prefer 3*N+2*N*NB)
893: // *
894: this._dgebrd.Run(N, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
895: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);
896: NWORK = IE + N;
897: // *
898: // * Perform bidiagonal SVD, computing singular values only
899: // * (Workspace: need N+BDSPAC)
900: // *
901: this._dbdsdc.Run("U", "N", N, ref S, offset_s, ref WORK, IE + o_work, ref DUM, offset_dum
902: , 1, ref DUM, offset_dum, 1, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work
903: , ref IWORK, offset_iwork, ref INFO);
904: // *
905: }
906: else
907: {
908: if (WNTQO)
909: {
910: // *
911: // * Path 2 (M much larger than N, JOBZ = 'O')
912: // * N left singular vectors to be overwritten on A and
913: // * N right singular vectors to be computed in VT
914: // *
915: IR = 1;
916: // *
917: // * WORK(IR) is LDWRKR by N
918: // *
919: if (LWORK >= LDA * N + N * N + 3 * N + BDSPAC)
920: {
921: LDWRKR = LDA;
922: }
923: else
924: {
925: LDWRKR = (LWORK - N * N - 3 * N - BDSPAC) / N;
926: }
927: ITAU = IR + LDWRKR * N;
928: NWORK = ITAU + N;
929: // *
930: // * Compute A=Q*R
931: // * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
932: // *
933: this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, NWORK + o_work
934: , LWORK - NWORK + 1, ref IERR);
935: // *
936: // * Copy R to WORK(IR), zeroing out below it
937: // *
938: this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref WORK, IR + o_work
939: , LDWRKR);
940: this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref WORK, IR + 1 + o_work
941: , LDWRKR);
942: // *
943: // * Generate Q in A
944: // * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
945: // *
946: this._dorgqr.Run(M, N, N, ref A, offset_a, LDA, WORK, ITAU + o_work
947: , ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);
948: IE = ITAU;
949: ITAUQ = IE + N;
950: ITAUP = ITAUQ + N;
951: NWORK = ITAUP + N;
952: // *
953: // * Bidiagonalize R in VT, copying result to WORK(IR)
954: // * (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
955: // *
956: this._dgebrd.Run(N, N, ref WORK, IR + o_work, LDWRKR, ref S, offset_s, ref WORK, IE + o_work
957: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);
958: // *
959: // * WORK(IU) is N by N
960: // *
961: IU = NWORK;
962: NWORK = IU + N * N;
963: // *
964: // * Perform bidiagonal SVD, computing left singular vectors
965: // * of bidiagonal matrix in WORK(IU) and computing right
966: // * singular vectors of bidiagonal matrix in VT
967: // * (Workspace: need N+N*N+BDSPAC)
968: // *
969: this._dbdsdc.Run("U", "I", N, ref S, offset_s, ref WORK, IE + o_work, ref WORK, IU + o_work
970: , N, ref VT, offset_vt, LDVT, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work
971: , ref IWORK, offset_iwork, ref INFO);
972: // *
973: // * Overwrite WORK(IU) by left singular vectors of R
974: // * and VT by right singular vectors of R
975: // * (Workspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB)
976: // *
977: this._dormbr.Run("Q", "L", "N", N, N, N
978: , ref WORK, IR + o_work, LDWRKR, WORK, ITAUQ + o_work, ref WORK, IU + o_work, N, ref WORK, NWORK + o_work
979: , LWORK - NWORK + 1, ref IERR);
980: this._dormbr.Run("P", "R", "T", N, N, N
981: , ref WORK, IR + o_work, LDWRKR, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, NWORK + o_work
982: , LWORK - NWORK + 1, ref IERR);
983: // *
984: // * Multiply Q in A by left singular vectors of R in
985: // * WORK(IU), storing result in WORK(IR) and copying to A
986: // * (Workspace: need 2*N*N, prefer N*N+M*N)
987: // *
988: for (I = 1; (LDWRKR >= 0) ? (I <= M) : (I >= M); I += LDWRKR)
989: {
990: CHUNK = Math.Min(M - I + 1, LDWRKR);
991: this._dgemm.Run("N", "N", CHUNK, N, N, ONE
992: , A, I+1 * LDA + o_a, LDA, WORK, IU + o_work, N, ZERO, ref WORK, IR + o_work
993: , LDWRKR);
994: this._dlacpy.Run("F", CHUNK, N, WORK, IR + o_work, LDWRKR, ref A, I+1 * LDA + o_a
995: , LDA);
996: }
997: // *
998: }
999: else
1000: {
1001: if (WNTQS)
1002: {
1003: // *
1004: // * Path 3 (M much larger than N, JOBZ='S')
1005: // * N left singular vectors to be computed in U and
1006: // * N right singular vectors to be computed in VT
1007: // *
1008: IR = 1;
1009: // *
1010: // * WORK(IR) is N by N
1011: // *
1012: LDWRKR = N;
1013: ITAU = IR + LDWRKR * N;
1014: NWORK = ITAU + N;
1015: // *
1016: // * Compute A=Q*R
1017: // * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1018: // *
1019: this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, NWORK + o_work
1020: , LWORK - NWORK + 1, ref IERR);
1021: // *
1022: // * Copy R to WORK(IR), zeroing out below it
1023: // *
1024: this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref WORK, IR + o_work
1025: , LDWRKR);
1026: this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref WORK, IR + 1 + o_work
1027: , LDWRKR);
1028: // *
1029: // * Generate Q in A
1030: // * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1031: // *
1032: this._dorgqr.Run(M, N, N, ref A, offset_a, LDA, WORK, ITAU + o_work
1033: , ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);
1034: IE = ITAU;
1035: ITAUQ = IE + N;
1036: ITAUP = ITAUQ + N;
1037: NWORK = ITAUP + N;
1038: // *
1039: // * Bidiagonalize R in WORK(IR)
1040: // * (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
1041: // *
1042: this._dgebrd.Run(N, N, ref WORK, IR + o_work, LDWRKR, ref S, offset_s, ref WORK, IE + o_work
1043: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);
1044: // *
1045: // * Perform bidiagonal SVD, computing left singular vectors
1046: // * of bidiagoal matrix in U and computing right singular
1047: // * vectors of bidiagonal matrix in VT
1048: // * (Workspace: need N+BDSPAC)
1049: // *
1050: this._dbdsdc.Run("U", "I", N, ref S, offset_s, ref WORK, IE + o_work, ref U, offset_u
1051: , LDU, ref VT, offset_vt, LDVT, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work
1052: , ref IWORK, offset_iwork, ref INFO);
1053: // *
1054: // * Overwrite U by left singular vectors of R and VT
1055: // * by right singular vectors of R
1056: // * (Workspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1057: // *
1058: this._dormbr.Run("Q", "L", "N", N, N, N
1059: , ref WORK, IR + o_work, LDWRKR, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, NWORK + o_work
1060: , LWORK - NWORK + 1, ref IERR);
1061: // *
1062: this._dormbr.Run("P", "R", "T", N, N, N
1063: , ref WORK, IR + o_work, LDWRKR, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, NWORK + o_work
1064: , LWORK - NWORK + 1, ref IERR);
1065: // *
1066: // * Multiply Q in A by left singular vectors of R in
1067: // * WORK(IR), storing result in U
1068: // * (Workspace: need N*N)
1069: // *
1070: this._dlacpy.Run("F", N, N, U, offset_u, LDU, ref WORK, IR + o_work
1071: , LDWRKR);
1072: this._dgemm.Run("N", "N", M, N, N, ONE
1073: , A, offset_a, LDA, WORK, IR + o_work, LDWRKR, ZERO, ref U, offset_u
1074: , LDU);
1075: // *
1076: }
1077: else
1078: {
1079: if (WNTQA)
1080: {
1081: // *
1082: // * Path 4 (M much larger than N, JOBZ='A')
1083: // * M left singular vectors to be computed in U and
1084: // * N right singular vectors to be computed in VT
1085: // *
1086: IU = 1;
1087: // *
1088: // * WORK(IU) is N by N
1089: // *
1090: LDWRKU = N;
1091: ITAU = IU + LDWRKU * N;
1092: NWORK = ITAU + N;
1093: // *
1094: // * Compute A=Q*R, copying result to U
1095: // * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1096: // *
1097: this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, NWORK + o_work
1098: , LWORK - NWORK + 1, ref IERR);
1099: this._dlacpy.Run("L", M, N, A, offset_a, LDA, ref U, offset_u
1100: , LDU);
1101: // *
1102: // * Generate Q in U
1103: // * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1104: this._dorgqr.Run(M, M, N, ref U, offset_u, LDU, WORK, ITAU + o_work
1105: , ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);
1106: // *
1107: // * Produce R in A, zeroing out other entries
1108: // *
1109: this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref A, 2+1 * LDA + o_a
1110: , LDA);
1111: IE = ITAU;
1112: ITAUQ = IE + N;
1113: ITAUP = ITAUQ + N;
1114: NWORK = ITAUP + N;
1115: // *
1116: // * Bidiagonalize R in A
1117: // * (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
1118: // *
1119: this._dgebrd.Run(N, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
1120: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);
1121: // *
1122: // * Perform bidiagonal SVD, computing left singular vectors
1123: // * of bidiagonal matrix in WORK(IU) and computing right
1124: // * singular vectors of bidiagonal matrix in VT
1125: // * (Workspace: need N+N*N+BDSPAC)
1126: // *
1127: this._dbdsdc.Run("U", "I", N, ref S, offset_s, ref WORK, IE + o_work, ref WORK, IU + o_work
1128: , N, ref VT, offset_vt, LDVT, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work
1129: , ref IWORK, offset_iwork, ref INFO);
1130: // *
1131: // * Overwrite WORK(IU) by left singular vectors of R and VT
1132: // * by right singular vectors of R
1133: // * (Workspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1134: // *
1135: this._dormbr.Run("Q", "L", "N", N, N, N
1136: , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref WORK, IU + o_work, LDWRKU, ref WORK, NWORK + o_work
1137: , LWORK - NWORK + 1, ref IERR);
1138: this._dormbr.Run("P", "R", "T", N, N, N
1139: , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, NWORK + o_work
1140: , LWORK - NWORK + 1, ref IERR);
1141: // *
1142: // * Multiply Q in U by left singular vectors of R in
1143: // * WORK(IU), storing result in A
1144: // * (Workspace: need N*N)
1145: // *
1146: this._dgemm.Run("N", "N", M, N, N, ONE
1147: , U, offset_u, LDU, WORK, IU + o_work, LDWRKU, ZERO, ref A, offset_a
1148: , LDA);
1149: // *
1150: // * Copy left singular vectors of A from A to U
1151: // *
1152: this._dlacpy.Run("F", M, N, A, offset_a, LDA, ref U, offset_u
1153: , LDU);
1154: // *
1155: }
1156: }
1157: }
1158: }
1159: // *
1160: }
1161: else
1162: {
1163: // *
1164: // * M .LT. MNTHR
1165: // *
1166: // * Path 5 (M at least N, but not much larger)
1167: // * Reduce to bidiagonal form without QR decomposition
1168: // *
1169: IE = 1;
1170: ITAUQ = IE + N;
1171: ITAUP = ITAUQ + N;
1172: NWORK = ITAUP + N;
1173: // *
1174: // * Bidiagonalize A
1175: // * (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB)
1176: // *
1177: this._dgebrd.Run(M, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
1178: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);
1179: if (WNTQN)
1180: {
1181: // *
1182: // * Perform bidiagonal SVD, only computing singular values
1183: // * (Workspace: need N+BDSPAC)
1184: // *
1185: this._dbdsdc.Run("U", "N", N, ref S, offset_s, ref WORK, IE + o_work, ref DUM, offset_dum
1186: , 1, ref DUM, offset_dum, 1, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work
1187: , ref IWORK, offset_iwork, ref INFO);
1188: }
1189: else
1190: {
1191: if (WNTQO)
1192: {
1193: IU = NWORK;
1194: if (LWORK >= M * N + 3 * N + BDSPAC)
1195: {
1196: // *
1197: // * WORK( IU ) is M by N
1198: // *
1199: LDWRKU = M;
1200: NWORK = IU + LDWRKU * N;
1201: this._dlaset.Run("F", M, N, ZERO, ZERO, ref WORK, IU + o_work
1202: , LDWRKU);
1203: }
1204: else
1205: {
1206: // *
1207: // * WORK( IU ) is N by N
1208: // *
1209: LDWRKU = N;
1210: NWORK = IU + LDWRKU * N;
1211: // *
1212: // * WORK(IR) is LDWRKR by N
1213: // *
1214: IR = NWORK;
1215: LDWRKR = (LWORK - N * N - 3 * N) / N;
1216: }
1217: NWORK = IU + LDWRKU * N;
1218: // *
1219: // * Perform bidiagonal SVD, computing left singular vectors
1220: // * of bidiagonal matrix in WORK(IU) and computing right
1221: // * singular vectors of bidiagonal matrix in VT
1222: // * (Workspace: need N+N*N+BDSPAC)
1223: // *
1224: this._dbdsdc.Run("U", "I", N, ref S, offset_s, ref WORK, IE + o_work, ref WORK, IU + o_work
1225: , LDWRKU, ref VT, offset_vt, LDVT, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work
1226: , ref IWORK, offset_iwork, ref INFO);
1227: // *
1228: // * Overwrite VT by right singular vectors of A
1229: // * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1230: // *
1231: this._dormbr.Run("P", "R", "T", N, N, N
1232: , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, NWORK + o_work
1233: , LWORK - NWORK + 1, ref IERR);
1234: // *
1235: if (LWORK >= M * N + 3 * N + BDSPAC)
1236: {
1237: // *
1238: // * Overwrite WORK(IU) by left singular vectors of A
1239: // * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1240: // *
1241: this._dormbr.Run("Q", "L", "N", M, N, N
1242: , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref WORK, IU + o_work, LDWRKU, ref WORK, NWORK + o_work
1243: , LWORK - NWORK + 1, ref IERR);
1244: // *
1245: // * Copy left singular vectors of A from WORK(IU) to A
1246: // *
1247: this._dlacpy.Run("F", M, N, WORK, IU + o_work, LDWRKU, ref A, offset_a
1248: , LDA);
1249: }
1250: else
1251: {
1252: // *
1253: // * Generate Q in A
1254: // * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1255: // *
1256: this._dorgbr.Run("Q", M, N, N, ref A, offset_a, LDA
1257: , WORK, ITAUQ + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);
1258: // *
1259: // * Multiply Q in A by left singular vectors of
1260: // * bidiagonal matrix in WORK(IU), storing result in
1261: // * WORK(IR) and copying to A
1262: // * (Workspace: need 2*N*N, prefer N*N+M*N)
1263: // *
1264: for (I = 1; (LDWRKR >= 0) ? (I <= M) : (I >= M); I += LDWRKR)
1265: {
1266: CHUNK = Math.Min(M - I + 1, LDWRKR);
1267: this._dgemm.Run("N", "N", CHUNK, N, N, ONE
1268: , A, I+1 * LDA + o_a, LDA, WORK, IU + o_work, LDWRKU, ZERO, ref WORK, IR + o_work
1269: , LDWRKR);
1270: this._dlacpy.Run("F", CHUNK, N, WORK, IR + o_work, LDWRKR, ref A, I+1 * LDA + o_a
1271: , LDA);
1272: }
1273: }
1274: // *
1275: }
1276: else
1277: {
1278: if (WNTQS)
1279: {
1280: // *
1281: // * Perform bidiagonal SVD, computing left singular vectors
1282: // * of bidiagonal matrix in U and computing right singular
1283: // * vectors of bidiagonal matrix in VT
1284: // * (Workspace: need N+BDSPAC)
1285: // *
1286: this._dlaset.Run("F", M, N, ZERO, ZERO, ref U, offset_u
1287: , LDU);
1288: this._dbdsdc.Run("U", "I", N, ref S, offset_s, ref WORK, IE + o_work, ref U, offset_u
1289: , LDU, ref VT, offset_vt, LDVT, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work
1290: , ref IWORK, offset_iwork, ref INFO);
1291: // *
1292: // * Overwrite U by left singular vectors of A and VT
1293: // * by right singular vectors of A
1294: // * (Workspace: need 3*N, prefer 2*N+N*NB)
1295: // *
1296: this._dormbr.Run("Q", "L", "N", M, N, N
1297: , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, NWORK + o_work
1298: , LWORK - NWORK + 1, ref IERR);
1299: this._dormbr.Run("P", "R", "T", N, N, N
1300: , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, NWORK + o_work
1301: , LWORK - NWORK + 1, ref IERR);
1302: }
1303: else
1304: {
1305: if (WNTQA)
1306: {
1307: // *
1308: // * Perform bidiagonal SVD, computing left singular vectors
1309: // * of bidiagonal matrix in U and computing right singular
1310: // * vectors of bidiagonal matrix in VT
1311: // * (Workspace: need N+BDSPAC)
1312: // *
1313: this._dlaset.Run("F", M, M, ZERO, ZERO, ref U, offset_u
1314: , LDU);
1315: this._dbdsdc.Run("U", "I", N, ref S, offset_s, ref WORK, IE + o_work, ref U, offset_u
1316: , LDU, ref VT, offset_vt, LDVT, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work
1317: , ref IWORK, offset_iwork, ref INFO);
1318: // *
1319: // * Set the right corner of U to identity matrix
1320: // *
1321: if (M > N)
1322: {
1323: this._dlaset.Run("F", M - N, M - N, ZERO, ONE, ref U, N + 1+(N + 1) * LDU + o_u
1324: , LDU);
1325: }
1326: // *
1327: // * Overwrite U by left singular vectors of A and VT
1328: // * by right singular vectors of A
1329: // * (Workspace: need N*N+2*N+M, prefer N*N+2*N+M*NB)
1330: // *
1331: this._dormbr.Run("Q", "L", "N", M, M, N
1332: , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, NWORK + o_work
1333: , LWORK - NWORK + 1, ref IERR);
1334: this._dormbr.Run("P", "R", "T", N, N, M
1335: , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, NWORK + o_work
1336: , LWORK - NWORK + 1, ref IERR);
1337: }
1338: }
1339: }
1340: }
1341: // *
1342: }
1343: // *
1344: }
1345: else
1346: {
1347: // *
1348: // * A has more columns than rows. If A has sufficiently more
1349: // * columns than rows, first reduce using the LQ decomposition (if
1350: // * sufficient workspace available)
1351: // *
1352: if (N >= MNTHR)
1353: {
1354: // *
1355: if (WNTQN)
1356: {
1357: // *
1358: // * Path 1t (N much larger than M, JOBZ='N')
1359: // * No singular vectors to be computed
1360: // *
1361: ITAU = 1;
1362: NWORK = ITAU + M;
1363: // *
1364: // * Compute A=L*Q
1365: // * (Workspace: need 2*M, prefer M+M*NB)
1366: // *
1367: this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, NWORK + o_work
1368: , LWORK - NWORK + 1, ref IERR);
1369: // *
1370: // * Zero out above L
1371: // *
1372: this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref A, 1+2 * LDA + o_a
1373: , LDA);
1374: IE = 1;
1375: ITAUQ = IE + M;
1376: ITAUP = ITAUQ + M;
1377: NWORK = ITAUP + M;
1378: // *
1379: // * Bidiagonalize L in A
1380: // * (Workspace: need 4*M, prefer 3*M+2*M*NB)
1381: // *
1382: this._dgebrd.Run(M, M, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
1383: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);
1384: NWORK = IE + M;
1385: // *
1386: // * Perform bidiagonal SVD, computing singular values only
1387: // * (Workspace: need M+BDSPAC)
1388: // *
1389: this._dbdsdc.Run("U", "N", M, ref S, offset_s, ref WORK, IE + o_work, ref DUM, offset_dum
1390: , 1, ref DUM, offset_dum, 1, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work
1391: , ref IWORK, offset_iwork, ref INFO);
1392: // *
1393: }
1394: else
1395: {
1396: if (WNTQO)
1397: {
1398: // *
1399: // * Path 2t (N much larger than M, JOBZ='O')
1400: // * M right singular vectors to be overwritten on A and
1401: // * M left singular vectors to be computed in U
1402: // *
1403: IVT = 1;
1404: // *
1405: // * IVT is M by M
1406: // *
1407: IL = IVT + M * M;
1408: if (LWORK >= M * N + M * M + 3 * M + BDSPAC)
1409: {
1410: // *
1411: // * WORK(IL) is M by N
1412: // *
1413: LDWRKL = M;
1414: CHUNK = N;
1415: }
1416: else
1417: {
1418: LDWRKL = M;
1419: CHUNK = (LWORK - M * M) / M;
1420: }
1421: ITAU = IL + LDWRKL * M;
1422: NWORK = ITAU + M;
1423: // *
1424: // * Compute A=L*Q
1425: // * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1426: // *
1427: this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, NWORK + o_work
1428: , LWORK - NWORK + 1, ref IERR);
1429: // *
1430: // * Copy L to WORK(IL), zeroing about above it
1431: // *
1432: this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref WORK, IL + o_work
1433: , LDWRKL);
1434: this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref WORK, IL + LDWRKL + o_work
1435: , LDWRKL);
1436: // *
1437: // * Generate Q in A
1438: // * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1439: // *
1440: this._dorglq.Run(M, N, M, ref A, offset_a, LDA, WORK, ITAU + o_work
1441: , ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);
1442: IE = ITAU;
1443: ITAUQ = IE + M;
1444: ITAUP = ITAUQ + M;
1445: NWORK = ITAUP + M;
1446: // *
1447: // * Bidiagonalize L in WORK(IL)
1448: // * (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
1449: // *
1450: this._dgebrd.Run(M, M, ref WORK, IL + o_work, LDWRKL, ref S, offset_s, ref WORK, IE + o_work
1451: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);
1452: // *
1453: // * Perform bidiagonal SVD, computing left singular vectors
1454: // * of bidiagonal matrix in U, and computing right singular
1455: // * vectors of bidiagonal matrix in WORK(IVT)
1456: // * (Workspace: need M+M*M+BDSPAC)
1457: // *
1458: this._dbdsdc.Run("U", "I", M, ref S, offset_s, ref WORK, IE + o_work, ref U, offset_u
1459: , LDU, ref WORK, IVT + o_work, M, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work
1460: , ref IWORK, offset_iwork, ref INFO);
1461: // *
1462: // * Overwrite U by left singular vectors of L and WORK(IVT)
1463: // * by right singular vectors of L
1464: // * (Workspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB)
1465: // *
1466: this._dormbr.Run("Q", "L", "N", M, M, M
1467: , ref WORK, IL + o_work, LDWRKL, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, NWORK + o_work
1468: , LWORK - NWORK + 1, ref IERR);
1469: this._dormbr.Run("P", "R", "T", M, M, M
1470: , ref WORK, IL + o_work, LDWRKL, WORK, ITAUP + o_work, ref WORK, IVT + o_work, M, ref WORK, NWORK + o_work
1471: , LWORK - NWORK + 1, ref IERR);
1472: // *
1473: // * Multiply right singular vectors of L in WORK(IVT) by Q
1474: // * in A, storing result in WORK(IL) and copying to A
1475: // * (Workspace: need 2*M*M, prefer M*M+M*N)
1476: // *
1477: for (I = 1; (CHUNK >= 0) ? (I <= N) : (I >= N); I += CHUNK)
1478: {
1479: BLK = Math.Min(N - I + 1, CHUNK);
1480: this._dgemm.Run("N", "N", M, BLK, M, ONE
1481: , WORK, IVT + o_work, M, A, 1+I * LDA + o_a, LDA, ZERO, ref WORK, IL + o_work
1482: , LDWRKL);
1483: this._dlacpy.Run("F", M, BLK, WORK, IL + o_work, LDWRKL, ref A, 1+I * LDA + o_a
1484: , LDA);
1485: }
1486: // *
1487: }
1488: else
1489: {
1490: if (WNTQS)
1491: {
1492: // *
1493: // * Path 3t (N much larger than M, JOBZ='S')
1494: // * M right singular vectors to be computed in VT and
1495: // * M left singular vectors to be computed in U
1496: // *
1497: IL = 1;
1498: // *
1499: // * WORK(IL) is M by M
1500: // *
1501: LDWRKL = M;
1502: ITAU = IL + LDWRKL * M;
1503: NWORK = ITAU + M;
1504: // *
1505: // * Compute A=L*Q
1506: // * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1507: // *
1508: this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, NWORK + o_work
1509: , LWORK - NWORK + 1, ref IERR);
1510: // *
1511: // * Copy L to WORK(IL), zeroing out above it
1512: // *
1513: this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref WORK, IL + o_work
1514: , LDWRKL);
1515: this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref WORK, IL + LDWRKL + o_work
1516: , LDWRKL);
1517: // *
1518: // * Generate Q in A
1519: // * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1520: // *
1521: this._dorglq.Run(M, N, M, ref A, offset_a, LDA, WORK, ITAU + o_work
1522: , ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);
1523: IE = ITAU;
1524: ITAUQ = IE + M;
1525: ITAUP = ITAUQ + M;
1526: NWORK = ITAUP + M;
1527: // *
1528: // * Bidiagonalize L in WORK(IU), copying result to U
1529: // * (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
1530: // *
1531: this._dgebrd.Run(M, M, ref WORK, IL + o_work, LDWRKL, ref S, offset_s, ref WORK, IE + o_work
1532: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);
1533: // *
1534: // * Perform bidiagonal SVD, computing left singular vectors
1535: // * of bidiagonal matrix in U and computing right singular
1536: // * vectors of bidiagonal matrix in VT
1537: // * (Workspace: need M+BDSPAC)
1538: // *
1539: this._dbdsdc.Run("U", "I", M, ref S, offset_s, ref WORK, IE + o_work, ref U, offset_u
1540: , LDU, ref VT, offset_vt, LDVT, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work
1541: , ref IWORK, offset_iwork, ref INFO);
1542: // *
1543: // * Overwrite U by left singular vectors of L and VT
1544: // * by right singular vectors of L
1545: // * (Workspace: need M*M+3*M, prefer M*M+2*M+M*NB)
1546: // *
1547: this._dormbr.Run("Q", "L", "N", M, M, M
1548: , ref WORK, IL + o_work, LDWRKL, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, NWORK + o_work
1549: , LWORK - NWORK + 1, ref IERR);
1550: this._dormbr.Run("P", "R", "T", M, M, M
1551: , ref WORK, IL + o_work, LDWRKL, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, NWORK + o_work
1552: , LWORK - NWORK + 1, ref IERR);
1553: // *
1554: // * Multiply right singular vectors of L in WORK(IL) by
1555: // * Q in A, storing result in VT
1556: // * (Workspace: need M*M)
1557: // *
1558: this._dlacpy.Run("F", M, M, VT, offset_vt, LDVT, ref WORK, IL + o_work
1559: , LDWRKL);
1560: this._dgemm.Run("N", "N", M, N, M, ONE
1561: , WORK, IL + o_work, LDWRKL, A, offset_a, LDA, ZERO, ref VT, offset_vt
1562: , LDVT);
1563: // *
1564: }
1565: else
1566: {
1567: if (WNTQA)
1568: {
1569: // *
1570: // * Path 4t (N much larger than M, JOBZ='A')
1571: // * N right singular vectors to be computed in VT and
1572: // * M left singular vectors to be computed in U
1573: // *
1574: IVT = 1;
1575: // *
1576: // * WORK(IVT) is M by M
1577: // *
1578: LDWKVT = M;
1579: ITAU = IVT + LDWKVT * M;
1580: NWORK = ITAU + M;
1581: // *
1582: // * Compute A=L*Q, copying result to VT
1583: // * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1584: // *
1585: this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, NWORK + o_work
1586: , LWORK - NWORK + 1, ref IERR);
1587: this._dlacpy.Run("U", M, N, A, offset_a, LDA, ref VT, offset_vt
1588: , LDVT);
1589: // *
1590: // * Generate Q in VT
1591: // * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1592: // *
1593: this._dorglq.Run(N, N, M, ref VT, offset_vt, LDVT, WORK, ITAU + o_work
1594: , ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);
1595: // *
1596: // * Produce L in A, zeroing out other entries
1597: // *
1598: this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref A, 1+2 * LDA + o_a
1599: , LDA);
1600: IE = ITAU;
1601: ITAUQ = IE + M;
1602: ITAUP = ITAUQ + M;
1603: NWORK = ITAUP + M;
1604: // *
1605: // * Bidiagonalize L in A
1606: // * (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
1607: // *
1608: this._dgebrd.Run(M, M, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
1609: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);
1610: // *
1611: // * Perform bidiagonal SVD, computing left singular vectors
1612: // * of bidiagonal matrix in U and computing right singular
1613: // * vectors of bidiagonal matrix in WORK(IVT)
1614: // * (Workspace: need M+M*M+BDSPAC)
1615: // *
1616: this._dbdsdc.Run("U", "I", M, ref S, offset_s, ref WORK, IE + o_work, ref U, offset_u
1617: , LDU, ref WORK, IVT + o_work, LDWKVT, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work
1618: , ref IWORK, offset_iwork, ref INFO);
1619: // *
1620: // * Overwrite U by left singular vectors of L and WORK(IVT)
1621: // * by right singular vectors of L
1622: // * (Workspace: need M*M+3*M, prefer M*M+2*M+M*NB)
1623: // *
1624: this._dormbr.Run("Q", "L", "N", M, M, M
1625: , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, NWORK + o_work
1626: , LWORK - NWORK + 1, ref IERR);
1627: this._dormbr.Run("P", "R", "T", M, M, M
1628: , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref WORK, IVT + o_work, LDWKVT, ref WORK, NWORK + o_work
1629: , LWORK - NWORK + 1, ref IERR);
1630: // *
1631: // * Multiply right singular vectors of L in WORK(IVT) by
1632: // * Q in VT, storing result in A
1633: // * (Workspace: need M*M)
1634: // *
1635: this._dgemm.Run("N", "N", M, N, M, ONE
1636: , WORK, IVT + o_work, LDWKVT, VT, offset_vt, LDVT, ZERO, ref A, offset_a
1637: , LDA);
1638: // *
1639: // * Copy right singular vectors of A from A to VT
1640: // *
1641: this._dlacpy.Run("F", M, N, A, offset_a, LDA, ref VT, offset_vt
1642: , LDVT);
1643: // *
1644: }
1645: }
1646: }
1647: }
1648: // *
1649: }
1650: else
1651: {
1652: // *
1653: // * N .LT. MNTHR
1654: // *
1655: // * Path 5t (N greater than M, but not much larger)
1656: // * Reduce to bidiagonal form without LQ decomposition
1657: // *
1658: IE = 1;
1659: ITAUQ = IE + M;
1660: ITAUP = ITAUQ + M;
1661: NWORK = ITAUP + M;
1662: // *
1663: // * Bidiagonalize A
1664: // * (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
1665: // *
1666: this._dgebrd.Run(M, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
1667: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);
1668: if (WNTQN)
1669: {
1670: // *
1671: // * Perform bidiagonal SVD, only computing singular values
1672: // * (Workspace: need M+BDSPAC)
1673: // *
1674: this._dbdsdc.Run("L", "N", M, ref S, offset_s, ref WORK, IE + o_work, ref DUM, offset_dum
1675: , 1, ref DUM, offset_dum, 1, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work
1676: , ref IWORK, offset_iwork, ref INFO);
1677: }
1678: else
1679: {
1680: if (WNTQO)
1681: {
1682: LDWKVT = M;
1683: IVT = NWORK;
1684: if (LWORK >= M * N + 3 * M + BDSPAC)
1685: {
1686: // *
1687: // * WORK( IVT ) is M by N
1688: // *
1689: this._dlaset.Run("F", M, N, ZERO, ZERO, ref WORK, IVT + o_work
1690: , LDWKVT);
1691: NWORK = IVT + LDWKVT * N;
1692: }
1693: else
1694: {
1695: // *
1696: // * WORK( IVT ) is M by M
1697: // *
1698: NWORK = IVT + LDWKVT * M;
1699: IL = NWORK;
1700: // *
1701: // * WORK(IL) is M by CHUNK
1702: // *
1703: CHUNK = (LWORK - M * M - 3 * M) / M;
1704: }
1705: // *
1706: // * Perform bidiagonal SVD, computing left singular vectors
1707: // * of bidiagonal matrix in U and computing right singular
1708: // * vectors of bidiagonal matrix in WORK(IVT)
1709: // * (Workspace: need M*M+BDSPAC)
1710: // *
1711: this._dbdsdc.Run("L", "I", M, ref S, offset_s, ref WORK, IE + o_work, ref U, offset_u
1712: , LDU, ref WORK, IVT + o_work, LDWKVT, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work
1713: , ref IWORK, offset_iwork, ref INFO);
1714: // *
1715: // * Overwrite U by left singular vectors of A
1716: // * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1717: // *
1718: this._dormbr.Run("Q", "L", "N", M, M, N
1719: , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, NWORK + o_work
1720: , LWORK - NWORK + 1, ref IERR);
1721: // *
1722: if (LWORK >= M * N + 3 * M + BDSPAC)
1723: {
1724: // *
1725: // * Overwrite WORK(IVT) by left singular vectors of A
1726: // * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1727: // *
1728: this._dormbr.Run("P", "R", "T", M, N, M
1729: , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref WORK, IVT + o_work, LDWKVT, ref WORK, NWORK + o_work
1730: , LWORK - NWORK + 1, ref IERR);
1731: // *
1732: // * Copy right singular vectors of A from WORK(IVT) to A
1733: // *
1734: this._dlacpy.Run("F", M, N, WORK, IVT + o_work, LDWKVT, ref A, offset_a
1735: , LDA);
1736: }
1737: else
1738: {
1739: // *
1740: // * Generate P**T in A
1741: // * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1742: // *
1743: this._dorgbr.Run("P", M, N, M, ref A, offset_a, LDA
1744: , WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);
1745: // *
1746: // * Multiply Q in A by right singular vectors of
1747: // * bidiagonal matrix in WORK(IVT), storing result in
1748: // * WORK(IL) and copying to A
1749: // * (Workspace: need 2*M*M, prefer M*M+M*N)
1750: // *
1751: for (I = 1; (CHUNK >= 0) ? (I <= N) : (I >= N); I += CHUNK)
1752: {
1753: BLK = Math.Min(N - I + 1, CHUNK);
1754: this._dgemm.Run("N", "N", M, BLK, M, ONE
1755: , WORK, IVT + o_work, LDWKVT, A, 1+I * LDA + o_a, LDA, ZERO, ref WORK, IL + o_work
1756: , M);
1757: this._dlacpy.Run("F", M, BLK, WORK, IL + o_work, M, ref A, 1+I * LDA + o_a
1758: , LDA);
1759: }
1760: }
1761: }
1762: else
1763: {
1764: if (WNTQS)
1765: {
1766: // *
1767: // * Perform bidiagonal SVD, computing left singular vectors
1768: // * of bidiagonal matrix in U and computing right singular
1769: // * vectors of bidiagonal matrix in VT
1770: // * (Workspace: need M+BDSPAC)
1771: // *
1772: this._dlaset.Run("F", M, N, ZERO, ZERO, ref VT, offset_vt
1773: , LDVT);
1774: this._dbdsdc.Run("L", "I", M, ref S, offset_s, ref WORK, IE + o_work, ref U, offset_u
1775: , LDU, ref VT, offset_vt, LDVT, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work
1776: , ref IWORK, offset_iwork, ref INFO);
1777: // *
1778: // * Overwrite U by left singular vectors of A and VT
1779: // * by right singular vectors of A
1780: // * (Workspace: need 3*M, prefer 2*M+M*NB)
1781: // *
1782: this._dormbr.Run("Q", "L", "N", M, M, N
1783: , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, NWORK + o_work
1784: , LWORK - NWORK + 1, ref IERR);
1785: this._dormbr.Run("P", "R", "T", M, N, M
1786: , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, NWORK + o_work
1787: , LWORK - NWORK + 1, ref IERR);
1788: }
1789: else
1790: {
1791: if (WNTQA)
1792: {
1793: // *
1794: // * Perform bidiagonal SVD, computing left singular vectors
1795: // * of bidiagonal matrix in U and computing right singular
1796: // * vectors of bidiagonal matrix in VT
1797: // * (Workspace: need M+BDSPAC)
1798: // *
1799: this._dlaset.Run("F", N, N, ZERO, ZERO, ref VT, offset_vt
1800: , LDVT);
1801: this._dbdsdc.Run("L", "I", M, ref S, offset_s, ref WORK, IE + o_work, ref U, offset_u
1802: , LDU, ref VT, offset_vt, LDVT, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work
1803: , ref IWORK, offset_iwork, ref INFO);
1804: // *
1805: // * Set the right corner of VT to identity matrix
1806: // *
1807: if (N > M)
1808: {
1809: this._dlaset.Run("F", N - M, N - M, ZERO, ONE, ref VT, M + 1+(M + 1) * LDVT + o_vt
1810: , LDVT);
1811: }
1812: // *
1813: // * Overwrite U by left singular vectors of A and VT
1814: // * by right singular vectors of A
1815: // * (Workspace: need 2*M+N, prefer 2*M+N*NB)
1816: // *
1817: this._dormbr.Run("Q", "L", "N", M, M, N
1818: , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, NWORK + o_work
1819: , LWORK - NWORK + 1, ref IERR);
1820: this._dormbr.Run("P", "R", "T", N, N, M
1821: , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, NWORK + o_work
1822: , LWORK - NWORK + 1, ref IERR);
1823: }
1824: }
1825: }
1826: }
1827: // *
1828: }
1829: // *
1830: }
1831: // *
1832: // * Undo scaling if necessary
1833: // *
1834: if (ISCL == 1)
1835: {
1836: if (ANRM > BIGNUM)
1837: {
1838: this._dlascl.Run("G", 0, 0, BIGNUM, ANRM, MINMN
1839: , 1, ref S, offset_s, MINMN, ref IERR);
1840: }
1841: if (ANRM < SMLNUM)
1842: {
1843: this._dlascl.Run("G", 0, 0, SMLNUM, ANRM, MINMN
1844: , 1, ref S, offset_s, MINMN, ref IERR);
1845: }
1846: }
1847: // *
1848: // * Return optimal workspace in WORK(1)
1849: // *
1850: WORK[1 + o_work] = MAXWRK;
1851: // *
1852: return;
1853: // *
1854: // * End of DGESDD
1855: // *
1856:
1857: #endregion
1858:
1859: }
1860: }
1861: }