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: /// DGESVD computes the singular value decomposition (SVD) of a real
27: /// M-by-N matrix A, optionally computing the left and/or right singular
28: /// vectors. The SVD is written
29: ///
30: /// A = U * SIGMA * transpose(V)
31: ///
32: /// where SIGMA is an M-by-N matrix which is zero except for its
33: /// min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
34: /// V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA
35: /// are the singular values of A; they are real and non-negative, and
36: /// are returned in descending order. The first min(m,n) columns of
37: /// U and V are the left and right singular vectors of A.
38: ///
39: /// Note that the routine returns V**T, not V.
40: ///
41: ///</summary>
42: public class DGESVD
43: {
44:
45:
46: #region Dependencies
47:
48: DBDSQR _dbdsqr; DGEBRD _dgebrd; DGELQF _dgelqf; DGEMM _dgemm; DGEQRF _dgeqrf; DLACPY _dlacpy; DLASCL _dlascl;
49: DLASET _dlaset;DORGBR _dorgbr; DORGLQ _dorglq; DORGQR _dorgqr; DORMBR _dormbr; XERBLA _xerbla; LSAME _lsame;
50: ILAENV _ilaenv;DLAMCH _dlamch; DLANGE _dlange;
51:
52: #endregion
53:
54:
55: #region Fields
56:
57: const double ZERO = 0.0E0; const double ONE = 1.0E0; bool LQUERY = false; bool WNTUA = false; bool WNTUAS = false;
58: bool WNTUN = false;bool WNTUO = false; bool WNTUS = false; bool WNTVA = false; bool WNTVAS = false; bool WNTVN = false;
59: bool WNTVO = false;bool WNTVS = false; int BDSPAC = 0; int BLK = 0; int CHUNK = 0; int I = 0; int IE = 0; int IERR = 0;
60: int IR = 0;int ISCL = 0; int ITAU = 0; int ITAUP = 0; int ITAUQ = 0; int IU = 0; int IWORK = 0; int LDWRKR = 0;
61: int LDWRKU = 0;int MAXWRK = 0; int MINMN = 0; int MINWRK = 0; int MNTHR = 0; int NCU = 0; int NCVT = 0; int NRU = 0;
62: int NRVT = 0;int WRKBL = 0; double ANRM = 0; double BIGNUM = 0; double EPS = 0; double SMLNUM = 0;
63: double[] DUM = new double[1]; int offset_dum = 0;
64:
65: #endregion
66:
67: public DGESVD(DBDSQR dbdsqr, DGEBRD dgebrd, DGELQF dgelqf, DGEMM dgemm, DGEQRF dgeqrf, DLACPY dlacpy, DLASCL dlascl, DLASET dlaset, DORGBR dorgbr, DORGLQ dorglq
68: , DORGQR dorgqr, DORMBR dormbr, XERBLA xerbla, LSAME lsame, ILAENV ilaenv, DLAMCH dlamch, DLANGE dlange)
69: {
70:
71:
72: #region Set Dependencies
73:
74: this._dbdsqr = dbdsqr; this._dgebrd = dgebrd; this._dgelqf = dgelqf; this._dgemm = dgemm; this._dgeqrf = dgeqrf;
75: this._dlacpy = dlacpy;this._dlascl = dlascl; this._dlaset = dlaset; this._dorgbr = dorgbr; this._dorglq = dorglq;
76: this._dorgqr = dorgqr;this._dormbr = dormbr; this._xerbla = xerbla; this._lsame = lsame; this._ilaenv = ilaenv;
77: this._dlamch = dlamch;this._dlange = dlange;
78:
79: #endregion
80:
81: }
82:
83: public DGESVD()
84: {
85:
86:
87: #region Dependencies (Initialization)
88:
89: LSAME lsame = new LSAME();
90: DLAMC3 dlamc3 = new DLAMC3();
91: DLAS2 dlas2 = new DLAS2();
92: DCOPY dcopy = new DCOPY();
93: XERBLA xerbla = new XERBLA();
94: DLASQ5 dlasq5 = new DLASQ5();
95: DLAZQ4 dlazq4 = new DLAZQ4();
96: IEEECK ieeeck = new IEEECK();
97: IPARMQ iparmq = new IPARMQ();
98: DROT drot = new DROT();
99: DSCAL dscal = new DSCAL();
100: DSWAP dswap = new DSWAP();
101: DLAPY2 dlapy2 = new DLAPY2();
102: DNRM2 dnrm2 = new DNRM2();
103: DLASSQ dlassq = new DLASSQ();
104: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
105: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
106: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
107: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
108: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
109: DLARTG dlartg = new DLARTG(dlamch);
110: DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
111: DLASQ6 dlasq6 = new DLASQ6(dlamch);
112: DLAZQ3 dlazq3 = new DLAZQ3(dlasq5, dlasq6, dlazq4, dlamch);
113: DLASRT dlasrt = new DLASRT(lsame, xerbla);
114: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
115: DLASQ2 dlasq2 = new DLASQ2(dlazq3, dlasrt, xerbla, dlamch, ilaenv);
116: DLASQ1 dlasq1 = new DLASQ1(dcopy, dlas2, dlascl, dlasq2, dlasrt, xerbla, dlamch);
117: DLASR dlasr = new DLASR(lsame, xerbla);
118: DLASV2 dlasv2 = new DLASV2(dlamch);
119: DBDSQR dbdsqr = new DBDSQR(lsame, dlamch, dlartg, dlas2, dlasq1, dlasr, dlasv2, drot, dscal, dswap
120: , xerbla);
121: DGEMV dgemv = new DGEMV(lsame, xerbla);
122: DGER dger = new DGER(xerbla);
123: DLARF dlarf = new DLARF(dgemv, dger, lsame);
124: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
125: DGEBD2 dgebd2 = new DGEBD2(dlarf, dlarfg, xerbla);
126: DGEMM dgemm = new DGEMM(lsame, xerbla);
127: DLABRD dlabrd = new DLABRD(dgemv, dlarfg, dscal);
128: DGEBRD dgebrd = new DGEBRD(dgebd2, dgemm, dlabrd, xerbla, ilaenv);
129: DGELQ2 dgelq2 = new DGELQ2(dlarf, dlarfg, xerbla);
130: DTRMM dtrmm = new DTRMM(lsame, xerbla);
131: DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
132: DTRMV dtrmv = new DTRMV(lsame, xerbla);
133: DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);
134: DGELQF dgelqf = new DGELQF(dgelq2, dlarfb, dlarft, xerbla, ilaenv);
135: DGEQR2 dgeqr2 = new DGEQR2(dlarf, dlarfg, xerbla);
136: DGEQRF dgeqrf = new DGEQRF(dgeqr2, dlarfb, dlarft, xerbla, ilaenv);
137: DLACPY dlacpy = new DLACPY(lsame);
138: DLASET dlaset = new DLASET(lsame);
139: DORGL2 dorgl2 = new DORGL2(dlarf, dscal, xerbla);
140: DORGLQ dorglq = new DORGLQ(dlarfb, dlarft, dorgl2, xerbla, ilaenv);
141: DORG2R dorg2r = new DORG2R(dlarf, dscal, xerbla);
142: DORGQR dorgqr = new DORGQR(dlarfb, dlarft, dorg2r, xerbla, ilaenv);
143: DORGBR dorgbr = new DORGBR(lsame, ilaenv, dorglq, dorgqr, xerbla);
144: DORML2 dorml2 = new DORML2(lsame, dlarf, xerbla);
145: DORMLQ dormlq = new DORMLQ(lsame, ilaenv, dlarfb, dlarft, dorml2, xerbla);
146: DORM2R dorm2r = new DORM2R(lsame, dlarf, xerbla);
147: DORMQR dormqr = new DORMQR(lsame, ilaenv, dlarfb, dlarft, dorm2r, xerbla);
148: DORMBR dormbr = new DORMBR(lsame, ilaenv, dormlq, dormqr, xerbla);
149: DLANGE dlange = new DLANGE(dlassq, lsame);
150:
151: #endregion
152:
153:
154: #region Set Dependencies
155:
156: this._dbdsqr = dbdsqr; this._dgebrd = dgebrd; this._dgelqf = dgelqf; this._dgemm = dgemm; this._dgeqrf = dgeqrf;
157: this._dlacpy = dlacpy;this._dlascl = dlascl; this._dlaset = dlaset; this._dorgbr = dorgbr; this._dorglq = dorglq;
158: this._dorgqr = dorgqr;this._dormbr = dormbr; this._xerbla = xerbla; this._lsame = lsame; this._ilaenv = ilaenv;
159: this._dlamch = dlamch;this._dlange = dlange;
160:
161: #endregion
162:
163: }
164: /// <summary>
165: /// Purpose
166: /// =======
167: ///
168: /// DGESVD computes the singular value decomposition (SVD) of a real
169: /// M-by-N matrix A, optionally computing the left and/or right singular
170: /// vectors. The SVD is written
171: ///
172: /// A = U * SIGMA * transpose(V)
173: ///
174: /// where SIGMA is an M-by-N matrix which is zero except for its
175: /// min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
176: /// V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA
177: /// are the singular values of A; they are real and non-negative, and
178: /// are returned in descending order. The first min(m,n) columns of
179: /// U and V are the left and right singular vectors of A.
180: ///
181: /// Note that the routine returns V**T, not V.
182: ///
183: ///</summary>
184: /// <param name="JOBU">
185: /// (input) CHARACTER*1
186: /// Specifies options for computing all or part of the matrix U:
187: /// = 'A': all M columns of U are returned in array U:
188: /// = 'S': the first min(m,n) columns of U (the left singular
189: /// vectors) are returned in the array U;
190: /// = 'O': the first min(m,n) columns of U (the left singular
191: /// vectors) are overwritten on the array A;
192: /// = 'N': no columns of U (no left singular vectors) are
193: /// computed.
194: ///</param>
195: /// <param name="JOBVT">
196: /// (input) CHARACTER*1
197: /// Specifies options for computing all or part of the matrix
198: /// V**T:
199: /// = 'A': all N rows of V**T are returned in the array VT;
200: /// = 'S': the first min(m,n) rows of V**T (the right singular
201: /// vectors) are returned in the array VT;
202: /// = 'O': the first min(m,n) rows of V**T (the right singular
203: /// vectors) are overwritten on the array A;
204: /// = 'N': no rows of V**T (no right singular vectors) are
205: /// computed.
206: ///
207: /// JOBVT and JOBU cannot both be 'O'.
208: ///</param>
209: /// <param name="M">
210: /// (input) INTEGER
211: /// The number of rows of the input matrix A. M .GE. 0.
212: ///</param>
213: /// <param name="N">
214: /// (input) INTEGER
215: /// The number of columns of the input matrix A. N .GE. 0.
216: ///</param>
217: /// <param name="A">
218: /// = U * SIGMA * transpose(V)
219: ///</param>
220: /// <param name="LDA">
221: /// (input) INTEGER
222: /// The leading dimension of the array A. LDA .GE. max(1,M).
223: ///</param>
224: /// <param name="S">
225: /// (output) DOUBLE PRECISION array, dimension (min(M,N))
226: /// The singular values of A, sorted so that S(i) .GE. S(i+1).
227: ///</param>
228: /// <param name="U">
229: /// (output) DOUBLE PRECISION array, dimension (LDU,UCOL)
230: /// (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
231: /// If JOBU = 'A', U contains the M-by-M orthogonal matrix U;
232: /// if JOBU = 'S', U contains the first min(m,n) columns of U
233: /// (the left singular vectors, stored columnwise);
234: /// if JOBU = 'N' or 'O', U is not referenced.
235: ///</param>
236: /// <param name="LDU">
237: /// (input) INTEGER
238: /// The leading dimension of the array U. LDU .GE. 1; if
239: /// JOBU = 'S' or 'A', LDU .GE. M.
240: ///</param>
241: /// <param name="VT">
242: /// (output) DOUBLE PRECISION array, dimension (LDVT,N)
243: /// If JOBVT = 'A', VT contains the N-by-N orthogonal matrix
244: /// V**T;
245: /// if JOBVT = 'S', VT contains the first min(m,n) rows of
246: /// V**T (the right singular vectors, stored rowwise);
247: /// if JOBVT = 'N' or 'O', VT is not referenced.
248: ///</param>
249: /// <param name="LDVT">
250: /// (input) INTEGER
251: /// The leading dimension of the array VT. LDVT .GE. 1; if
252: /// JOBVT = 'A', LDVT .GE. N; if JOBVT = 'S', LDVT .GE. min(M,N).
253: ///</param>
254: /// <param name="WORK">
255: /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
256: /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
257: /// if INFO .GT. 0, WORK(2:MIN(M,N)) contains the unconverged
258: /// superdiagonal elements of an upper bidiagonal matrix B
259: /// whose diagonal is in S (not necessarily sorted). B
260: /// satisfies A = U * B * VT, so it has the same singular values
261: /// as A, and singular vectors related by U and VT.
262: ///</param>
263: /// <param name="LWORK">
264: /// (input) INTEGER
265: /// The dimension of the array WORK.
266: /// LWORK .GE. MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N)).
267: /// For good performance, LWORK should generally be larger.
268: ///
269: /// If LWORK = -1, then a workspace query is assumed; the routine
270: /// only calculates the optimal size of the WORK array, returns
271: /// this value as the first entry of the WORK array, and no error
272: /// message related to LWORK is issued by XERBLA.
273: ///</param>
274: /// <param name="INFO">
275: /// (output) INTEGER
276: /// = 0: successful exit.
277: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
278: /// .GT. 0: if DBDSQR did not converge, INFO specifies how many
279: /// superdiagonals of an intermediate bidiagonal form B
280: /// did not converge to zero. See the description of WORK
281: /// above for details.
282: ///</param>
283: public void Run(string JOBU, string JOBVT, int M, int N, ref double[] A, int offset_a, int LDA
284: , ref double[] S, int offset_s, ref double[] U, int offset_u, int LDU, ref double[] VT, int offset_vt, int LDVT, ref double[] WORK, int offset_work
285: , int LWORK, ref int INFO)
286: {
287:
288: #region Array Index Correction
289:
290: int o_a = -1 - LDA + offset_a; int o_s = -1 + offset_s; int o_u = -1 - LDU + offset_u;
291: int o_vt = -1 - LDVT + offset_vt; int o_work = -1 + offset_work;
292:
293: #endregion
294:
295:
296: #region Strings
297:
298: JOBU = JOBU.Substring(0, 1); JOBVT = JOBVT.Substring(0, 1);
299:
300: #endregion
301:
302:
303: #region Prolog
304:
305: // *
306: // * -- LAPACK driver routine (version 3.1) --
307: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
308: // * November 2006
309: // *
310: // * .. Scalar Arguments ..
311: // * ..
312: // * .. Array Arguments ..
313: // * ..
314: // *
315: // * Purpose
316: // * =======
317: // *
318: // * DGESVD computes the singular value decomposition (SVD) of a real
319: // * M-by-N matrix A, optionally computing the left and/or right singular
320: // * vectors. The SVD is written
321: // *
322: // * A = U * SIGMA * transpose(V)
323: // *
324: // * where SIGMA is an M-by-N matrix which is zero except for its
325: // * min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
326: // * V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA
327: // * are the singular values of A; they are real and non-negative, and
328: // * are returned in descending order. The first min(m,n) columns of
329: // * U and V are the left and right singular vectors of A.
330: // *
331: // * Note that the routine returns V**T, not V.
332: // *
333: // * Arguments
334: // * =========
335: // *
336: // * JOBU (input) CHARACTER*1
337: // * Specifies options for computing all or part of the matrix U:
338: // * = 'A': all M columns of U are returned in array U:
339: // * = 'S': the first min(m,n) columns of U (the left singular
340: // * vectors) are returned in the array U;
341: // * = 'O': the first min(m,n) columns of U (the left singular
342: // * vectors) are overwritten on the array A;
343: // * = 'N': no columns of U (no left singular vectors) are
344: // * computed.
345: // *
346: // * JOBVT (input) CHARACTER*1
347: // * Specifies options for computing all or part of the matrix
348: // * V**T:
349: // * = 'A': all N rows of V**T are returned in the array VT;
350: // * = 'S': the first min(m,n) rows of V**T (the right singular
351: // * vectors) are returned in the array VT;
352: // * = 'O': the first min(m,n) rows of V**T (the right singular
353: // * vectors) are overwritten on the array A;
354: // * = 'N': no rows of V**T (no right singular vectors) are
355: // * computed.
356: // *
357: // * JOBVT and JOBU cannot both be 'O'.
358: // *
359: // * M (input) INTEGER
360: // * The number of rows of the input matrix A. M >= 0.
361: // *
362: // * N (input) INTEGER
363: // * The number of columns of the input matrix A. N >= 0.
364: // *
365: // * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
366: // * On entry, the M-by-N matrix A.
367: // * On exit,
368: // * if JOBU = 'O', A is overwritten with the first min(m,n)
369: // * columns of U (the left singular vectors,
370: // * stored columnwise);
371: // * if JOBVT = 'O', A is overwritten with the first min(m,n)
372: // * rows of V**T (the right singular vectors,
373: // * stored rowwise);
374: // * if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A
375: // * are destroyed.
376: // *
377: // * LDA (input) INTEGER
378: // * The leading dimension of the array A. LDA >= max(1,M).
379: // *
380: // * S (output) DOUBLE PRECISION array, dimension (min(M,N))
381: // * The singular values of A, sorted so that S(i) >= S(i+1).
382: // *
383: // * U (output) DOUBLE PRECISION array, dimension (LDU,UCOL)
384: // * (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
385: // * If JOBU = 'A', U contains the M-by-M orthogonal matrix U;
386: // * if JOBU = 'S', U contains the first min(m,n) columns of U
387: // * (the left singular vectors, stored columnwise);
388: // * if JOBU = 'N' or 'O', U is not referenced.
389: // *
390: // * LDU (input) INTEGER
391: // * The leading dimension of the array U. LDU >= 1; if
392: // * JOBU = 'S' or 'A', LDU >= M.
393: // *
394: // * VT (output) DOUBLE PRECISION array, dimension (LDVT,N)
395: // * If JOBVT = 'A', VT contains the N-by-N orthogonal matrix
396: // * V**T;
397: // * if JOBVT = 'S', VT contains the first min(m,n) rows of
398: // * V**T (the right singular vectors, stored rowwise);
399: // * if JOBVT = 'N' or 'O', VT is not referenced.
400: // *
401: // * LDVT (input) INTEGER
402: // * The leading dimension of the array VT. LDVT >= 1; if
403: // * JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N).
404: // *
405: // * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
406: // * On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
407: // * if INFO > 0, WORK(2:MIN(M,N)) contains the unconverged
408: // * superdiagonal elements of an upper bidiagonal matrix B
409: // * whose diagonal is in S (not necessarily sorted). B
410: // * satisfies A = U * B * VT, so it has the same singular values
411: // * as A, and singular vectors related by U and VT.
412: // *
413: // * LWORK (input) INTEGER
414: // * The dimension of the array WORK.
415: // * LWORK >= MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N)).
416: // * For good performance, LWORK should generally be larger.
417: // *
418: // * If LWORK = -1, then a workspace query is assumed; the routine
419: // * only calculates the optimal size of the WORK array, returns
420: // * this value as the first entry of the WORK array, and no error
421: // * message related to LWORK is issued by XERBLA.
422: // *
423: // * INFO (output) INTEGER
424: // * = 0: successful exit.
425: // * < 0: if INFO = -i, the i-th argument had an illegal value.
426: // * > 0: if DBDSQR did not converge, INFO specifies how many
427: // * superdiagonals of an intermediate bidiagonal form B
428: // * did not converge to zero. See the description of WORK
429: // * above for details.
430: // *
431: // * =====================================================================
432: // *
433: // * .. Parameters ..
434: // * ..
435: // * .. Local Scalars ..
436: // * ..
437: // * .. Local Arrays ..
438: // * ..
439: // * .. External Subroutines ..
440: // * ..
441: // * .. External Functions ..
442: // * ..
443: // * .. Intrinsic Functions ..
444: // INTRINSIC MAX, MIN, SQRT;
445: // * ..
446: // * .. Executable Statements ..
447: // *
448: // * Test the input arguments
449: // *
450:
451: #endregion
452:
453:
454: #region Body
455:
456: INFO = 0;
457: MINMN = Math.Min(M, N);
458: WNTUA = this._lsame.Run(JOBU, "A");
459: WNTUS = this._lsame.Run(JOBU, "S");
460: WNTUAS = WNTUA || WNTUS;
461: WNTUO = this._lsame.Run(JOBU, "O");
462: WNTUN = this._lsame.Run(JOBU, "N");
463: WNTVA = this._lsame.Run(JOBVT, "A");
464: WNTVS = this._lsame.Run(JOBVT, "S");
465: WNTVAS = WNTVA || WNTVS;
466: WNTVO = this._lsame.Run(JOBVT, "O");
467: WNTVN = this._lsame.Run(JOBVT, "N");
468: LQUERY = (LWORK == - 1);
469: // *
470: if (!(WNTUA || WNTUS || WNTUO || WNTUN))
471: {
472: INFO = - 1;
473: }
474: else
475: {
476: if (!(WNTVA || WNTVS || WNTVO || WNTVN) || (WNTVO && WNTUO))
477: {
478: INFO = - 2;
479: }
480: else
481: {
482: if (M < 0)
483: {
484: INFO = - 3;
485: }
486: else
487: {
488: if (N < 0)
489: {
490: INFO = - 4;
491: }
492: else
493: {
494: if (LDA < Math.Max(1, M))
495: {
496: INFO = - 6;
497: }
498: else
499: {
500: if (LDU < 1 || (WNTUAS && LDU < M))
501: {
502: INFO = - 9;
503: }
504: else
505: {
506: if (LDVT < 1 || (WNTVA && LDVT < N) || (WNTVS && LDVT < MINMN))
507: {
508: INFO = - 11;
509: }
510: }
511: }
512: }
513: }
514: }
515: }
516: // *
517: // * Compute workspace
518: // * (Note: Comments in the code beginning "Workspace:" describe the
519: // * minimal amount of workspace needed at that point in the code,
520: // * as well as the preferred amount for good performance.
521: // * NB refers to the optimal block size for the immediately
522: // * following subroutine, as returned by ILAENV.)
523: // *
524: if (INFO == 0)
525: {
526: MINWRK = 1;
527: MAXWRK = 1;
528: if (M >= N && MINMN > 0)
529: {
530: // *
531: // * Compute space needed for DBDSQR
532: // *
533: MNTHR = this._ilaenv.Run(6, "DGESVD", JOBU + JOBVT, M, N, 0, 0);
534: BDSPAC = 5 * N;
535: if (M >= MNTHR)
536: {
537: if (WNTUN)
538: {
539: // *
540: // * Path 1 (M much larger than N, JOBU='N')
541: // *
542: MAXWRK = N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N, - 1, - 1);
543: MAXWRK = Math.Max(MAXWRK, 3 * N + 2 * N * this._ilaenv.Run(1, "DGEBRD", " ", N, N, - 1, - 1));
544: if (WNTVO || WNTVAS) MAXWRK = Math.Max(MAXWRK, 3 * N + (N - 1) * this._ilaenv.Run(1, "DORGBR", "P", N, N, N, - 1));
545: MAXWRK = Math.Max(MAXWRK, BDSPAC);
546: MINWRK = Math.Max(4 * N, BDSPAC);
547: }
548: else
549: {
550: if (WNTUO && WNTVN)
551: {
552: // *
553: // * Path 2 (M much larger than N, JOBU='O', JOBVT='N')
554: // *
555: WRKBL = N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N, - 1, - 1);
556: WRKBL = Math.Max(WRKBL, N + N * this._ilaenv.Run(1, "DORGQR", " ", M, N, N, - 1));
557: WRKBL = Math.Max(WRKBL, 3 * N + 2 * N * this._ilaenv.Run(1, "DGEBRD", " ", N, N, - 1, - 1));
558: WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORGBR", "Q", N, N, N, - 1));
559: WRKBL = Math.Max(WRKBL, BDSPAC);
560: MAXWRK = Math.Max(N * N + WRKBL, N * N + M * N + N);
561: MINWRK = Math.Max(3 * N + M, BDSPAC);
562: }
563: else
564: {
565: if (WNTUO && WNTVAS)
566: {
567: // *
568: // * Path 3 (M much larger than N, JOBU='O', JOBVT='S' or
569: // * 'A')
570: // *
571: WRKBL = N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N, - 1, - 1);
572: WRKBL = Math.Max(WRKBL, N + N * this._ilaenv.Run(1, "DORGQR", " ", M, N, N, - 1));
573: WRKBL = Math.Max(WRKBL, 3 * N + 2 * N * this._ilaenv.Run(1, "DGEBRD", " ", N, N, - 1, - 1));
574: WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORGBR", "Q", N, N, N, - 1));
575: WRKBL = Math.Max(WRKBL, 3 * N + (N - 1) * this._ilaenv.Run(1, "DORGBR", "P", N, N, N, - 1));
576: WRKBL = Math.Max(WRKBL, BDSPAC);
577: MAXWRK = Math.Max(N * N + WRKBL, N * N + M * N + N);
578: MINWRK = Math.Max(3 * N + M, BDSPAC);
579: }
580: else
581: {
582: if (WNTUS && WNTVN)
583: {
584: // *
585: // * Path 4 (M much larger than N, JOBU='S', JOBVT='N')
586: // *
587: WRKBL = N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N, - 1, - 1);
588: WRKBL = Math.Max(WRKBL, N + N * this._ilaenv.Run(1, "DORGQR", " ", M, N, N, - 1));
589: WRKBL = Math.Max(WRKBL, 3 * N + 2 * N * this._ilaenv.Run(1, "DGEBRD", " ", N, N, - 1, - 1));
590: WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORGBR", "Q", N, N, N, - 1));
591: WRKBL = Math.Max(WRKBL, BDSPAC);
592: MAXWRK = N * N + WRKBL;
593: MINWRK = Math.Max(3 * N + M, BDSPAC);
594: }
595: else
596: {
597: if (WNTUS && WNTVO)
598: {
599: // *
600: // * Path 5 (M much larger than N, JOBU='S', JOBVT='O')
601: // *
602: WRKBL = N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N, - 1, - 1);
603: WRKBL = Math.Max(WRKBL, N + N * this._ilaenv.Run(1, "DORGQR", " ", M, N, N, - 1));
604: WRKBL = Math.Max(WRKBL, 3 * N + 2 * N * this._ilaenv.Run(1, "DGEBRD", " ", N, N, - 1, - 1));
605: WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORGBR", "Q", N, N, N, - 1));
606: WRKBL = Math.Max(WRKBL, 3 * N + (N - 1) * this._ilaenv.Run(1, "DORGBR", "P", N, N, N, - 1));
607: WRKBL = Math.Max(WRKBL, BDSPAC);
608: MAXWRK = 2 * N * N + WRKBL;
609: MINWRK = Math.Max(3 * N + M, BDSPAC);
610: }
611: else
612: {
613: if (WNTUS && WNTVAS)
614: {
615: // *
616: // * Path 6 (M much larger than N, JOBU='S', JOBVT='S' or
617: // * 'A')
618: // *
619: WRKBL = N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N, - 1, - 1);
620: WRKBL = Math.Max(WRKBL, N + N * this._ilaenv.Run(1, "DORGQR", " ", M, N, N, - 1));
621: WRKBL = Math.Max(WRKBL, 3 * N + 2 * N * this._ilaenv.Run(1, "DGEBRD", " ", N, N, - 1, - 1));
622: WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORGBR", "Q", N, N, N, - 1));
623: WRKBL = Math.Max(WRKBL, 3 * N + (N - 1) * this._ilaenv.Run(1, "DORGBR", "P", N, N, N, - 1));
624: WRKBL = Math.Max(WRKBL, BDSPAC);
625: MAXWRK = N * N + WRKBL;
626: MINWRK = Math.Max(3 * N + M, BDSPAC);
627: }
628: else
629: {
630: if (WNTUA && WNTVN)
631: {
632: // *
633: // * Path 7 (M much larger than N, JOBU='A', JOBVT='N')
634: // *
635: WRKBL = N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N, - 1, - 1);
636: WRKBL = Math.Max(WRKBL, N + M * this._ilaenv.Run(1, "DORGQR", " ", M, M, N, - 1));
637: WRKBL = Math.Max(WRKBL, 3 * N + 2 * N * this._ilaenv.Run(1, "DGEBRD", " ", N, N, - 1, - 1));
638: WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORGBR", "Q", N, N, N, - 1));
639: WRKBL = Math.Max(WRKBL, BDSPAC);
640: MAXWRK = N * N + WRKBL;
641: MINWRK = Math.Max(3 * N + M, BDSPAC);
642: }
643: else
644: {
645: if (WNTUA && WNTVO)
646: {
647: // *
648: // * Path 8 (M much larger than N, JOBU='A', JOBVT='O')
649: // *
650: WRKBL = N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N, - 1, - 1);
651: WRKBL = Math.Max(WRKBL, N + M * this._ilaenv.Run(1, "DORGQR", " ", M, M, N, - 1));
652: WRKBL = Math.Max(WRKBL, 3 * N + 2 * N * this._ilaenv.Run(1, "DGEBRD", " ", N, N, - 1, - 1));
653: WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORGBR", "Q", N, N, N, - 1));
654: WRKBL = Math.Max(WRKBL, 3 * N + (N - 1) * this._ilaenv.Run(1, "DORGBR", "P", N, N, N, - 1));
655: WRKBL = Math.Max(WRKBL, BDSPAC);
656: MAXWRK = 2 * N * N + WRKBL;
657: MINWRK = Math.Max(3 * N + M, BDSPAC);
658: }
659: else
660: {
661: if (WNTUA && WNTVAS)
662: {
663: // *
664: // * Path 9 (M much larger than N, JOBU='A', JOBVT='S' or
665: // * 'A')
666: // *
667: WRKBL = N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N, - 1, - 1);
668: WRKBL = Math.Max(WRKBL, N + M * this._ilaenv.Run(1, "DORGQR", " ", M, M, N, - 1));
669: WRKBL = Math.Max(WRKBL, 3 * N + 2 * N * this._ilaenv.Run(1, "DGEBRD", " ", N, N, - 1, - 1));
670: WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORGBR", "Q", N, N, N, - 1));
671: WRKBL = Math.Max(WRKBL, 3 * N + (N - 1) * this._ilaenv.Run(1, "DORGBR", "P", N, N, N, - 1));
672: WRKBL = Math.Max(WRKBL, BDSPAC);
673: MAXWRK = N * N + WRKBL;
674: MINWRK = Math.Max(3 * N + M, BDSPAC);
675: }
676: }
677: }
678: }
679: }
680: }
681: }
682: }
683: }
684: }
685: else
686: {
687: // *
688: // * Path 10 (M at least N, but not much larger)
689: // *
690: MAXWRK = 3 * N + (M + N) * this._ilaenv.Run(1, "DGEBRD", " ", M, N, - 1, - 1);
691: if (WNTUS || WNTUO) MAXWRK = Math.Max(MAXWRK, 3 * N + N * this._ilaenv.Run(1, "DORGBR", "Q", M, N, N, - 1));
692: if (WNTUA) MAXWRK = Math.Max(MAXWRK, 3 * N + M * this._ilaenv.Run(1, "DORGBR", "Q", M, M, N, - 1));
693: if (!WNTVN) MAXWRK = Math.Max(MAXWRK, 3 * N + (N - 1) * this._ilaenv.Run(1, "DORGBR", "P", N, N, N, - 1));
694: MAXWRK = Math.Max(MAXWRK, BDSPAC);
695: MINWRK = Math.Max(3 * N + M, BDSPAC);
696: }
697: }
698: else
699: {
700: if (MINMN > 0)
701: {
702: // *
703: // * Compute space needed for DBDSQR
704: // *
705: MNTHR = this._ilaenv.Run(6, "DGESVD", JOBU + JOBVT, M, N, 0, 0);
706: BDSPAC = 5 * M;
707: if (N >= MNTHR)
708: {
709: if (WNTVN)
710: {
711: // *
712: // * Path 1t(N much larger than M, JOBVT='N')
713: // *
714: MAXWRK = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N, - 1, - 1);
715: MAXWRK = Math.Max(MAXWRK, 3 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M, - 1, - 1));
716: if (WNTUO || WNTUAS) MAXWRK = Math.Max(MAXWRK, 3 * M + M * this._ilaenv.Run(1, "DORGBR", "Q", M, M, M, - 1));
717: MAXWRK = Math.Max(MAXWRK, BDSPAC);
718: MINWRK = Math.Max(4 * M, BDSPAC);
719: }
720: else
721: {
722: if (WNTVO && WNTUN)
723: {
724: // *
725: // * Path 2t(N much larger than M, JOBU='N', JOBVT='O')
726: // *
727: WRKBL = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N, - 1, - 1);
728: WRKBL = Math.Max(WRKBL, M + M * this._ilaenv.Run(1, "DORGLQ", " ", M, N, M, - 1));
729: WRKBL = Math.Max(WRKBL, 3 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M, - 1, - 1));
730: WRKBL = Math.Max(WRKBL, 3 * M + (M - 1) * this._ilaenv.Run(1, "DORGBR", "P", M, M, M, - 1));
731: WRKBL = Math.Max(WRKBL, BDSPAC);
732: MAXWRK = Math.Max(M * M + WRKBL, M * M + M * N + M);
733: MINWRK = Math.Max(3 * M + N, BDSPAC);
734: }
735: else
736: {
737: if (WNTVO && WNTUAS)
738: {
739: // *
740: // * Path 3t(N much larger than M, JOBU='S' or 'A',
741: // * JOBVT='O')
742: // *
743: WRKBL = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N, - 1, - 1);
744: WRKBL = Math.Max(WRKBL, M + M * this._ilaenv.Run(1, "DORGLQ", " ", M, N, M, - 1));
745: WRKBL = Math.Max(WRKBL, 3 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M, - 1, - 1));
746: WRKBL = Math.Max(WRKBL, 3 * M + (M - 1) * this._ilaenv.Run(1, "DORGBR", "P", M, M, M, - 1));
747: WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORGBR", "Q", M, M, M, - 1));
748: WRKBL = Math.Max(WRKBL, BDSPAC);
749: MAXWRK = Math.Max(M * M + WRKBL, M * M + M * N + M);
750: MINWRK = Math.Max(3 * M + N, BDSPAC);
751: }
752: else
753: {
754: if (WNTVS && WNTUN)
755: {
756: // *
757: // * Path 4t(N much larger than M, JOBU='N', JOBVT='S')
758: // *
759: WRKBL = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N, - 1, - 1);
760: WRKBL = Math.Max(WRKBL, M + M * this._ilaenv.Run(1, "DORGLQ", " ", M, N, M, - 1));
761: WRKBL = Math.Max(WRKBL, 3 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M, - 1, - 1));
762: WRKBL = Math.Max(WRKBL, 3 * M + (M - 1) * this._ilaenv.Run(1, "DORGBR", "P", M, M, M, - 1));
763: WRKBL = Math.Max(WRKBL, BDSPAC);
764: MAXWRK = M * M + WRKBL;
765: MINWRK = Math.Max(3 * M + N, BDSPAC);
766: }
767: else
768: {
769: if (WNTVS && WNTUO)
770: {
771: // *
772: // * Path 5t(N much larger than M, JOBU='O', JOBVT='S')
773: // *
774: WRKBL = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N, - 1, - 1);
775: WRKBL = Math.Max(WRKBL, M + M * this._ilaenv.Run(1, "DORGLQ", " ", M, N, M, - 1));
776: WRKBL = Math.Max(WRKBL, 3 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M, - 1, - 1));
777: WRKBL = Math.Max(WRKBL, 3 * M + (M - 1) * this._ilaenv.Run(1, "DORGBR", "P", M, M, M, - 1));
778: WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORGBR", "Q", M, M, M, - 1));
779: WRKBL = Math.Max(WRKBL, BDSPAC);
780: MAXWRK = 2 * M * M + WRKBL;
781: MINWRK = Math.Max(3 * M + N, BDSPAC);
782: }
783: else
784: {
785: if (WNTVS && WNTUAS)
786: {
787: // *
788: // * Path 6t(N much larger than M, JOBU='S' or 'A',
789: // * JOBVT='S')
790: // *
791: WRKBL = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N, - 1, - 1);
792: WRKBL = Math.Max(WRKBL, M + M * this._ilaenv.Run(1, "DORGLQ", " ", M, N, M, - 1));
793: WRKBL = Math.Max(WRKBL, 3 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M, - 1, - 1));
794: WRKBL = Math.Max(WRKBL, 3 * M + (M - 1) * this._ilaenv.Run(1, "DORGBR", "P", M, M, M, - 1));
795: WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORGBR", "Q", M, M, M, - 1));
796: WRKBL = Math.Max(WRKBL, BDSPAC);
797: MAXWRK = M * M + WRKBL;
798: MINWRK = Math.Max(3 * M + N, BDSPAC);
799: }
800: else
801: {
802: if (WNTVA && WNTUN)
803: {
804: // *
805: // * Path 7t(N much larger than M, JOBU='N', JOBVT='A')
806: // *
807: WRKBL = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N, - 1, - 1);
808: WRKBL = Math.Max(WRKBL, M + N * this._ilaenv.Run(1, "DORGLQ", " ", N, N, M, - 1));
809: WRKBL = Math.Max(WRKBL, 3 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M, - 1, - 1));
810: WRKBL = Math.Max(WRKBL, 3 * M + (M - 1) * this._ilaenv.Run(1, "DORGBR", "P", M, M, M, - 1));
811: WRKBL = Math.Max(WRKBL, BDSPAC);
812: MAXWRK = M * M + WRKBL;
813: MINWRK = Math.Max(3 * M + N, BDSPAC);
814: }
815: else
816: {
817: if (WNTVA && WNTUO)
818: {
819: // *
820: // * Path 8t(N much larger than M, JOBU='O', JOBVT='A')
821: // *
822: WRKBL = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N, - 1, - 1);
823: WRKBL = Math.Max(WRKBL, M + N * this._ilaenv.Run(1, "DORGLQ", " ", N, N, M, - 1));
824: WRKBL = Math.Max(WRKBL, 3 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M, - 1, - 1));
825: WRKBL = Math.Max(WRKBL, 3 * M + (M - 1) * this._ilaenv.Run(1, "DORGBR", "P", M, M, M, - 1));
826: WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORGBR", "Q", M, M, M, - 1));
827: WRKBL = Math.Max(WRKBL, BDSPAC);
828: MAXWRK = 2 * M * M + WRKBL;
829: MINWRK = Math.Max(3 * M + N, BDSPAC);
830: }
831: else
832: {
833: if (WNTVA && WNTUAS)
834: {
835: // *
836: // * Path 9t(N much larger than M, JOBU='S' or 'A',
837: // * JOBVT='A')
838: // *
839: WRKBL = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N, - 1, - 1);
840: WRKBL = Math.Max(WRKBL, M + N * this._ilaenv.Run(1, "DORGLQ", " ", N, N, M, - 1));
841: WRKBL = Math.Max(WRKBL, 3 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M, - 1, - 1));
842: WRKBL = Math.Max(WRKBL, 3 * M + (M - 1) * this._ilaenv.Run(1, "DORGBR", "P", M, M, M, - 1));
843: WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORGBR", "Q", M, M, M, - 1));
844: WRKBL = Math.Max(WRKBL, BDSPAC);
845: MAXWRK = M * M + WRKBL;
846: MINWRK = Math.Max(3 * M + N, BDSPAC);
847: }
848: }
849: }
850: }
851: }
852: }
853: }
854: }
855: }
856: }
857: else
858: {
859: // *
860: // * Path 10t(N greater than M, but not much larger)
861: // *
862: MAXWRK = 3 * M + (M + N) * this._ilaenv.Run(1, "DGEBRD", " ", M, N, - 1, - 1);
863: if (WNTVS || WNTVO) MAXWRK = Math.Max(MAXWRK, 3 * M + M * this._ilaenv.Run(1, "DORGBR", "P", M, N, M, - 1));
864: if (WNTVA) MAXWRK = Math.Max(MAXWRK, 3 * M + N * this._ilaenv.Run(1, "DORGBR", "P", N, N, M, - 1));
865: if (!WNTUN) MAXWRK = Math.Max(MAXWRK, 3 * M + (M - 1) * this._ilaenv.Run(1, "DORGBR", "Q", M, M, M, - 1));
866: MAXWRK = Math.Max(MAXWRK, BDSPAC);
867: MINWRK = Math.Max(3 * M + N, BDSPAC);
868: }
869: }
870: }
871: MAXWRK = Math.Max(MAXWRK, MINWRK);
872: WORK[1 + o_work] = MAXWRK;
873: // *
874: if (LWORK < MINWRK && !LQUERY)
875: {
876: INFO = - 13;
877: }
878: }
879: // *
880: if (INFO != 0)
881: {
882: this._xerbla.Run("DGESVD", - INFO);
883: return;
884: }
885: else
886: {
887: if (LQUERY)
888: {
889: return;
890: }
891: }
892: // *
893: // * Quick return if possible
894: // *
895: if (M == 0 || N == 0)
896: {
897: return;
898: }
899: // *
900: // * Get machine constants
901: // *
902: EPS = this._dlamch.Run("P");
903: SMLNUM = Math.Sqrt(this._dlamch.Run("S")) / EPS;
904: BIGNUM = ONE / SMLNUM;
905: // *
906: // * Scale A if max element outside range [SMLNUM,BIGNUM]
907: // *
908: ANRM = this._dlange.Run("M", M, N, A, offset_a, LDA, ref DUM, offset_dum);
909: ISCL = 0;
910: if (ANRM > ZERO && ANRM < SMLNUM)
911: {
912: ISCL = 1;
913: this._dlascl.Run("G", 0, 0, ANRM, SMLNUM, M
914: , N, ref A, offset_a, LDA, ref IERR);
915: }
916: else
917: {
918: if (ANRM > BIGNUM)
919: {
920: ISCL = 1;
921: this._dlascl.Run("G", 0, 0, ANRM, BIGNUM, M
922: , N, ref A, offset_a, LDA, ref IERR);
923: }
924: }
925: // *
926: if (M >= N)
927: {
928: // *
929: // * A has at least as many rows as columns. If A has sufficiently
930: // * more rows than columns, first reduce using the QR
931: // * decomposition (if sufficient workspace available)
932: // *
933: if (M >= MNTHR)
934: {
935: // *
936: if (WNTUN)
937: {
938: // *
939: // * Path 1 (M much larger than N, JOBU='N')
940: // * No left singular vectors to be computed
941: // *
942: ITAU = 1;
943: IWORK = ITAU + N;
944: // *
945: // * Compute A=Q*R
946: // * (Workspace: need 2*N, prefer N+N*NB)
947: // *
948: this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
949: , LWORK - IWORK + 1, ref IERR);
950: // *
951: // * Zero out below R
952: // *
953: this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref A, 2+1 * LDA + o_a
954: , LDA);
955: IE = 1;
956: ITAUQ = IE + N;
957: ITAUP = ITAUQ + N;
958: IWORK = ITAUP + N;
959: // *
960: // * Bidiagonalize R in A
961: // * (Workspace: need 4*N, prefer 3*N+2*N*NB)
962: // *
963: this._dgebrd.Run(N, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
964: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
965: NCVT = 0;
966: if (WNTVO || WNTVAS)
967: {
968: // *
969: // * If right singular vectors desired, generate P'.
970: // * (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
971: // *
972: this._dorgbr.Run("P", N, N, N, ref A, offset_a, LDA
973: , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
974: NCVT = N;
975: }
976: IWORK = IE + N;
977: // *
978: // * Perform bidiagonal QR iteration, computing right
979: // * singular vectors of A in A if desired
980: // * (Workspace: need BDSPAC)
981: // *
982: this._dbdsqr.Run("U", N, NCVT, 0, 0, ref S, offset_s
983: , ref WORK, IE + o_work, ref A, offset_a, LDA, ref DUM, offset_dum, 1, ref DUM, offset_dum
984: , 1, ref WORK, IWORK + o_work, ref INFO);
985: // *
986: // * If right singular vectors desired in VT, copy them there
987: // *
988: if (WNTVAS)
989: {
990: this._dlacpy.Run("F", N, N, A, offset_a, LDA, ref VT, offset_vt
991: , LDVT);
992: }
993: // *
994: }
995: else
996: {
997: if (WNTUO && WNTVN)
998: {
999: // *
1000: // * Path 2 (M much larger than N, JOBU='O', JOBVT='N')
1001: // * N left singular vectors to be overwritten on A and
1002: // * no right singular vectors to be computed
1003: // *
1004: if (LWORK >= N * N + Math.Max(4 * N, BDSPAC))
1005: {
1006: // *
1007: // * Sufficient workspace for a fast algorithm
1008: // *
1009: IR = 1;
1010: if (LWORK >= Math.Max(WRKBL, LDA * N + N) + LDA * N)
1011: {
1012: // *
1013: // * WORK(IU) is LDA by N, WORK(IR) is LDA by N
1014: // *
1015: LDWRKU = LDA;
1016: LDWRKR = LDA;
1017: }
1018: else
1019: {
1020: if (LWORK >= Math.Max(WRKBL, LDA * N + N) + N * N)
1021: {
1022: // *
1023: // * WORK(IU) is LDA by N, WORK(IR) is N by N
1024: // *
1025: LDWRKU = LDA;
1026: LDWRKR = N;
1027: }
1028: else
1029: {
1030: // *
1031: // * WORK(IU) is LDWRKU by N, WORK(IR) is N by N
1032: // *
1033: LDWRKU = (LWORK - N * N - N) / N;
1034: LDWRKR = N;
1035: }
1036: }
1037: ITAU = IR + LDWRKR * N;
1038: IWORK = ITAU + N;
1039: // *
1040: // * Compute A=Q*R
1041: // * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1042: // *
1043: this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
1044: , LWORK - IWORK + 1, ref IERR);
1045: // *
1046: // * Copy R to WORK(IR) and zero out below it
1047: // *
1048: this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref WORK, IR + o_work
1049: , LDWRKR);
1050: this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref WORK, IR + 1 + o_work
1051: , LDWRKR);
1052: // *
1053: // * Generate Q in A
1054: // * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1055: // *
1056: this._dorgqr.Run(M, N, N, ref A, offset_a, LDA, WORK, ITAU + o_work
1057: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1058: IE = ITAU;
1059: ITAUQ = IE + N;
1060: ITAUP = ITAUQ + N;
1061: IWORK = ITAUP + N;
1062: // *
1063: // * Bidiagonalize R in WORK(IR)
1064: // * (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
1065: // *
1066: this._dgebrd.Run(N, N, ref WORK, IR + o_work, LDWRKR, ref S, offset_s, ref WORK, IE + o_work
1067: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1068: // *
1069: // * Generate left vectors bidiagonalizing R
1070: // * (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
1071: // *
1072: this._dorgbr.Run("Q", N, N, N, ref WORK, IR + o_work, LDWRKR
1073: , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1074: IWORK = IE + N;
1075: // *
1076: // * Perform bidiagonal QR iteration, computing left
1077: // * singular vectors of R in WORK(IR)
1078: // * (Workspace: need N*N+BDSPAC)
1079: // *
1080: this._dbdsqr.Run("U", N, 0, N, 0, ref S, offset_s
1081: , ref WORK, IE + o_work, ref DUM, offset_dum, 1, ref WORK, IR + o_work, LDWRKR, ref DUM, offset_dum
1082: , 1, ref WORK, IWORK + o_work, ref INFO);
1083: IU = IE + N;
1084: // *
1085: // * Multiply Q in A by left singular vectors of R in
1086: // * WORK(IR), storing result in WORK(IU) and copying to A
1087: // * (Workspace: need N*N+2*N, prefer N*N+M*N+N)
1088: // *
1089: for (I = 1; (LDWRKU >= 0) ? (I <= M) : (I >= M); I += LDWRKU)
1090: {
1091: CHUNK = Math.Min(M - I + 1, LDWRKU);
1092: this._dgemm.Run("N", "N", CHUNK, N, N, ONE
1093: , A, I+1 * LDA + o_a, LDA, WORK, IR + o_work, LDWRKR, ZERO, ref WORK, IU + o_work
1094: , LDWRKU);
1095: this._dlacpy.Run("F", CHUNK, N, WORK, IU + o_work, LDWRKU, ref A, I+1 * LDA + o_a
1096: , LDA);
1097: }
1098: // *
1099: }
1100: else
1101: {
1102: // *
1103: // * Insufficient workspace for a fast algorithm
1104: // *
1105: IE = 1;
1106: ITAUQ = IE + N;
1107: ITAUP = ITAUQ + N;
1108: IWORK = ITAUP + N;
1109: // *
1110: // * Bidiagonalize A
1111: // * (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB)
1112: // *
1113: this._dgebrd.Run(M, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
1114: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1115: // *
1116: // * Generate left vectors bidiagonalizing A
1117: // * (Workspace: need 4*N, prefer 3*N+N*NB)
1118: // *
1119: this._dorgbr.Run("Q", M, N, N, ref A, offset_a, LDA
1120: , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1121: IWORK = IE + N;
1122: // *
1123: // * Perform bidiagonal QR iteration, computing left
1124: // * singular vectors of A in A
1125: // * (Workspace: need BDSPAC)
1126: // *
1127: this._dbdsqr.Run("U", N, 0, M, 0, ref S, offset_s
1128: , ref WORK, IE + o_work, ref DUM, offset_dum, 1, ref A, offset_a, LDA, ref DUM, offset_dum
1129: , 1, ref WORK, IWORK + o_work, ref INFO);
1130: // *
1131: }
1132: // *
1133: }
1134: else
1135: {
1136: if (WNTUO && WNTVAS)
1137: {
1138: // *
1139: // * Path 3 (M much larger than N, JOBU='O', JOBVT='S' or 'A')
1140: // * N left singular vectors to be overwritten on A and
1141: // * N right singular vectors to be computed in VT
1142: // *
1143: if (LWORK >= N * N + Math.Max(4 * N, BDSPAC))
1144: {
1145: // *
1146: // * Sufficient workspace for a fast algorithm
1147: // *
1148: IR = 1;
1149: if (LWORK >= Math.Max(WRKBL, LDA * N + N) + LDA * N)
1150: {
1151: // *
1152: // * WORK(IU) is LDA by N and WORK(IR) is LDA by N
1153: // *
1154: LDWRKU = LDA;
1155: LDWRKR = LDA;
1156: }
1157: else
1158: {
1159: if (LWORK >= Math.Max(WRKBL, LDA * N + N) + N * N)
1160: {
1161: // *
1162: // * WORK(IU) is LDA by N and WORK(IR) is N by N
1163: // *
1164: LDWRKU = LDA;
1165: LDWRKR = N;
1166: }
1167: else
1168: {
1169: // *
1170: // * WORK(IU) is LDWRKU by N and WORK(IR) is N by N
1171: // *
1172: LDWRKU = (LWORK - N * N - N) / N;
1173: LDWRKR = N;
1174: }
1175: }
1176: ITAU = IR + LDWRKR * N;
1177: IWORK = ITAU + N;
1178: // *
1179: // * Compute A=Q*R
1180: // * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1181: // *
1182: this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
1183: , LWORK - IWORK + 1, ref IERR);
1184: // *
1185: // * Copy R to VT, zeroing out below it
1186: // *
1187: this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref VT, offset_vt
1188: , LDVT);
1189: if (N > 1)
1190: {
1191: this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref VT, 2+1 * LDVT + o_vt
1192: , LDVT);
1193: }
1194: // *
1195: // * Generate Q in A
1196: // * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1197: // *
1198: this._dorgqr.Run(M, N, N, ref A, offset_a, LDA, WORK, ITAU + o_work
1199: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1200: IE = ITAU;
1201: ITAUQ = IE + N;
1202: ITAUP = ITAUQ + N;
1203: IWORK = ITAUP + N;
1204: // *
1205: // * Bidiagonalize R in VT, copying result to WORK(IR)
1206: // * (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
1207: // *
1208: this._dgebrd.Run(N, N, ref VT, offset_vt, LDVT, ref S, offset_s, ref WORK, IE + o_work
1209: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1210: this._dlacpy.Run("L", N, N, VT, offset_vt, LDVT, ref WORK, IR + o_work
1211: , LDWRKR);
1212: // *
1213: // * Generate left vectors bidiagonalizing R in WORK(IR)
1214: // * (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
1215: // *
1216: this._dorgbr.Run("Q", N, N, N, ref WORK, IR + o_work, LDWRKR
1217: , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1218: // *
1219: // * Generate right vectors bidiagonalizing R in VT
1220: // * (Workspace: need N*N+4*N-1, prefer N*N+3*N+(N-1)*NB)
1221: // *
1222: this._dorgbr.Run("P", N, N, N, ref VT, offset_vt, LDVT
1223: , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1224: IWORK = IE + N;
1225: // *
1226: // * Perform bidiagonal QR iteration, computing left
1227: // * singular vectors of R in WORK(IR) and computing right
1228: // * singular vectors of R in VT
1229: // * (Workspace: need N*N+BDSPAC)
1230: // *
1231: this._dbdsqr.Run("U", N, N, N, 0, ref S, offset_s
1232: , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref WORK, IR + o_work, LDWRKR, ref DUM, offset_dum
1233: , 1, ref WORK, IWORK + o_work, ref INFO);
1234: IU = IE + N;
1235: // *
1236: // * Multiply Q in A by left singular vectors of R in
1237: // * WORK(IR), storing result in WORK(IU) and copying to A
1238: // * (Workspace: need N*N+2*N, prefer N*N+M*N+N)
1239: // *
1240: for (I = 1; (LDWRKU >= 0) ? (I <= M) : (I >= M); I += LDWRKU)
1241: {
1242: CHUNK = Math.Min(M - I + 1, LDWRKU);
1243: this._dgemm.Run("N", "N", CHUNK, N, N, ONE
1244: , A, I+1 * LDA + o_a, LDA, WORK, IR + o_work, LDWRKR, ZERO, ref WORK, IU + o_work
1245: , LDWRKU);
1246: this._dlacpy.Run("F", CHUNK, N, WORK, IU + o_work, LDWRKU, ref A, I+1 * LDA + o_a
1247: , LDA);
1248: }
1249: // *
1250: }
1251: else
1252: {
1253: // *
1254: // * Insufficient workspace for a fast algorithm
1255: // *
1256: ITAU = 1;
1257: IWORK = ITAU + N;
1258: // *
1259: // * Compute A=Q*R
1260: // * (Workspace: need 2*N, prefer N+N*NB)
1261: // *
1262: this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
1263: , LWORK - IWORK + 1, ref IERR);
1264: // *
1265: // * Copy R to VT, zeroing out below it
1266: // *
1267: this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref VT, offset_vt
1268: , LDVT);
1269: if (N > 1)
1270: {
1271: this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref VT, 2+1 * LDVT + o_vt
1272: , LDVT);
1273: }
1274: // *
1275: // * Generate Q in A
1276: // * (Workspace: need 2*N, prefer N+N*NB)
1277: // *
1278: this._dorgqr.Run(M, N, N, ref A, offset_a, LDA, WORK, ITAU + o_work
1279: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1280: IE = ITAU;
1281: ITAUQ = IE + N;
1282: ITAUP = ITAUQ + N;
1283: IWORK = ITAUP + N;
1284: // *
1285: // * Bidiagonalize R in VT
1286: // * (Workspace: need 4*N, prefer 3*N+2*N*NB)
1287: // *
1288: this._dgebrd.Run(N, N, ref VT, offset_vt, LDVT, ref S, offset_s, ref WORK, IE + o_work
1289: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1290: // *
1291: // * Multiply Q in A by left vectors bidiagonalizing R
1292: // * (Workspace: need 3*N+M, prefer 3*N+M*NB)
1293: // *
1294: this._dormbr.Run("Q", "R", "N", M, N, N
1295: , ref VT, offset_vt, LDVT, WORK, ITAUQ + o_work, ref A, offset_a, LDA, ref WORK, IWORK + o_work
1296: , LWORK - IWORK + 1, ref IERR);
1297: // *
1298: // * Generate right vectors bidiagonalizing R in VT
1299: // * (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
1300: // *
1301: this._dorgbr.Run("P", N, N, N, ref VT, offset_vt, LDVT
1302: , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1303: IWORK = IE + N;
1304: // *
1305: // * Perform bidiagonal QR iteration, computing left
1306: // * singular vectors of A in A and computing right
1307: // * singular vectors of A in VT
1308: // * (Workspace: need BDSPAC)
1309: // *
1310: this._dbdsqr.Run("U", N, N, M, 0, ref S, offset_s
1311: , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref A, offset_a, LDA, ref DUM, offset_dum
1312: , 1, ref WORK, IWORK + o_work, ref INFO);
1313: // *
1314: }
1315: // *
1316: }
1317: else
1318: {
1319: if (WNTUS)
1320: {
1321: // *
1322: if (WNTVN)
1323: {
1324: // *
1325: // * Path 4 (M much larger than N, JOBU='S', JOBVT='N')
1326: // * N left singular vectors to be computed in U and
1327: // * no right singular vectors to be computed
1328: // *
1329: if (LWORK >= N * N + Math.Max(4 * N, BDSPAC))
1330: {
1331: // *
1332: // * Sufficient workspace for a fast algorithm
1333: // *
1334: IR = 1;
1335: if (LWORK >= WRKBL + LDA * N)
1336: {
1337: // *
1338: // * WORK(IR) is LDA by N
1339: // *
1340: LDWRKR = LDA;
1341: }
1342: else
1343: {
1344: // *
1345: // * WORK(IR) is N by N
1346: // *
1347: LDWRKR = N;
1348: }
1349: ITAU = IR + LDWRKR * N;
1350: IWORK = ITAU + N;
1351: // *
1352: // * Compute A=Q*R
1353: // * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1354: // *
1355: this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
1356: , LWORK - IWORK + 1, ref IERR);
1357: // *
1358: // * Copy R to WORK(IR), zeroing out below it
1359: // *
1360: this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref WORK, IR + o_work
1361: , LDWRKR);
1362: this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref WORK, IR + 1 + o_work
1363: , LDWRKR);
1364: // *
1365: // * Generate Q in A
1366: // * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1367: // *
1368: this._dorgqr.Run(M, N, N, ref A, offset_a, LDA, WORK, ITAU + o_work
1369: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1370: IE = ITAU;
1371: ITAUQ = IE + N;
1372: ITAUP = ITAUQ + N;
1373: IWORK = ITAUP + N;
1374: // *
1375: // * Bidiagonalize R in WORK(IR)
1376: // * (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
1377: // *
1378: this._dgebrd.Run(N, N, ref WORK, IR + o_work, LDWRKR, ref S, offset_s, ref WORK, IE + o_work
1379: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1380: // *
1381: // * Generate left vectors bidiagonalizing R in WORK(IR)
1382: // * (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
1383: // *
1384: this._dorgbr.Run("Q", N, N, N, ref WORK, IR + o_work, LDWRKR
1385: , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1386: IWORK = IE + N;
1387: // *
1388: // * Perform bidiagonal QR iteration, computing left
1389: // * singular vectors of R in WORK(IR)
1390: // * (Workspace: need N*N+BDSPAC)
1391: // *
1392: this._dbdsqr.Run("U", N, 0, N, 0, ref S, offset_s
1393: , ref WORK, IE + o_work, ref DUM, offset_dum, 1, ref WORK, IR + o_work, LDWRKR, ref DUM, offset_dum
1394: , 1, ref WORK, IWORK + o_work, ref INFO);
1395: // *
1396: // * Multiply Q in A by left singular vectors of R in
1397: // * WORK(IR), storing result in U
1398: // * (Workspace: need N*N)
1399: // *
1400: this._dgemm.Run("N", "N", M, N, N, ONE
1401: , A, offset_a, LDA, WORK, IR + o_work, LDWRKR, ZERO, ref U, offset_u
1402: , LDU);
1403: // *
1404: }
1405: else
1406: {
1407: // *
1408: // * Insufficient workspace for a fast algorithm
1409: // *
1410: ITAU = 1;
1411: IWORK = ITAU + N;
1412: // *
1413: // * Compute A=Q*R, copying result to U
1414: // * (Workspace: need 2*N, prefer N+N*NB)
1415: // *
1416: this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
1417: , LWORK - IWORK + 1, ref IERR);
1418: this._dlacpy.Run("L", M, N, A, offset_a, LDA, ref U, offset_u
1419: , LDU);
1420: // *
1421: // * Generate Q in U
1422: // * (Workspace: need 2*N, prefer N+N*NB)
1423: // *
1424: this._dorgqr.Run(M, N, N, ref U, offset_u, LDU, WORK, ITAU + o_work
1425: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1426: IE = ITAU;
1427: ITAUQ = IE + N;
1428: ITAUP = ITAUQ + N;
1429: IWORK = ITAUP + N;
1430: // *
1431: // * Zero out below R in A
1432: // *
1433: this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref A, 2+1 * LDA + o_a
1434: , LDA);
1435: // *
1436: // * Bidiagonalize R in A
1437: // * (Workspace: need 4*N, prefer 3*N+2*N*NB)
1438: // *
1439: this._dgebrd.Run(N, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
1440: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1441: // *
1442: // * Multiply Q in U by left vectors bidiagonalizing R
1443: // * (Workspace: need 3*N+M, prefer 3*N+M*NB)
1444: // *
1445: this._dormbr.Run("Q", "R", "N", M, N, N
1446: , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, IWORK + o_work
1447: , LWORK - IWORK + 1, ref IERR);
1448: IWORK = IE + N;
1449: // *
1450: // * Perform bidiagonal QR iteration, computing left
1451: // * singular vectors of A in U
1452: // * (Workspace: need BDSPAC)
1453: // *
1454: this._dbdsqr.Run("U", N, 0, M, 0, ref S, offset_s
1455: , ref WORK, IE + o_work, ref DUM, offset_dum, 1, ref U, offset_u, LDU, ref DUM, offset_dum
1456: , 1, ref WORK, IWORK + o_work, ref INFO);
1457: // *
1458: }
1459: // *
1460: }
1461: else
1462: {
1463: if (WNTVO)
1464: {
1465: // *
1466: // * Path 5 (M much larger than N, JOBU='S', JOBVT='O')
1467: // * N left singular vectors to be computed in U and
1468: // * N right singular vectors to be overwritten on A
1469: // *
1470: if (LWORK >= 2 * N * N + Math.Max(4 * N, BDSPAC))
1471: {
1472: // *
1473: // * Sufficient workspace for a fast algorithm
1474: // *
1475: IU = 1;
1476: if (LWORK >= WRKBL + 2 * LDA * N)
1477: {
1478: // *
1479: // * WORK(IU) is LDA by N and WORK(IR) is LDA by N
1480: // *
1481: LDWRKU = LDA;
1482: IR = IU + LDWRKU * N;
1483: LDWRKR = LDA;
1484: }
1485: else
1486: {
1487: if (LWORK >= WRKBL + (LDA + N) * N)
1488: {
1489: // *
1490: // * WORK(IU) is LDA by N and WORK(IR) is N by N
1491: // *
1492: LDWRKU = LDA;
1493: IR = IU + LDWRKU * N;
1494: LDWRKR = N;
1495: }
1496: else
1497: {
1498: // *
1499: // * WORK(IU) is N by N and WORK(IR) is N by N
1500: // *
1501: LDWRKU = N;
1502: IR = IU + LDWRKU * N;
1503: LDWRKR = N;
1504: }
1505: }
1506: ITAU = IR + LDWRKR * N;
1507: IWORK = ITAU + N;
1508: // *
1509: // * Compute A=Q*R
1510: // * (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
1511: // *
1512: this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
1513: , LWORK - IWORK + 1, ref IERR);
1514: // *
1515: // * Copy R to WORK(IU), zeroing out below it
1516: // *
1517: this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref WORK, IU + o_work
1518: , LDWRKU);
1519: this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref WORK, IU + 1 + o_work
1520: , LDWRKU);
1521: // *
1522: // * Generate Q in A
1523: // * (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
1524: // *
1525: this._dorgqr.Run(M, N, N, ref A, offset_a, LDA, WORK, ITAU + o_work
1526: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1527: IE = ITAU;
1528: ITAUQ = IE + N;
1529: ITAUP = ITAUQ + N;
1530: IWORK = ITAUP + N;
1531: // *
1532: // * Bidiagonalize R in WORK(IU), copying result to
1533: // * WORK(IR)
1534: // * (Workspace: need 2*N*N+4*N,
1535: // * prefer 2*N*N+3*N+2*N*NB)
1536: // *
1537: this._dgebrd.Run(N, N, ref WORK, IU + o_work, LDWRKU, ref S, offset_s, ref WORK, IE + o_work
1538: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1539: this._dlacpy.Run("U", N, N, WORK, IU + o_work, LDWRKU, ref WORK, IR + o_work
1540: , LDWRKR);
1541: // *
1542: // * Generate left bidiagonalizing vectors in WORK(IU)
1543: // * (Workspace: need 2*N*N+4*N, prefer 2*N*N+3*N+N*NB)
1544: // *
1545: this._dorgbr.Run("Q", N, N, N, ref WORK, IU + o_work, LDWRKU
1546: , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1547: // *
1548: // * Generate right bidiagonalizing vectors in WORK(IR)
1549: // * (Workspace: need 2*N*N+4*N-1,
1550: // * prefer 2*N*N+3*N+(N-1)*NB)
1551: // *
1552: this._dorgbr.Run("P", N, N, N, ref WORK, IR + o_work, LDWRKR
1553: , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1554: IWORK = IE + N;
1555: // *
1556: // * Perform bidiagonal QR iteration, computing left
1557: // * singular vectors of R in WORK(IU) and computing
1558: // * right singular vectors of R in WORK(IR)
1559: // * (Workspace: need 2*N*N+BDSPAC)
1560: // *
1561: this._dbdsqr.Run("U", N, N, N, 0, ref S, offset_s
1562: , ref WORK, IE + o_work, ref WORK, IR + o_work, LDWRKR, ref WORK, IU + o_work, LDWRKU, ref DUM, offset_dum
1563: , 1, ref WORK, IWORK + o_work, ref INFO);
1564: // *
1565: // * Multiply Q in A by left singular vectors of R in
1566: // * WORK(IU), storing result in U
1567: // * (Workspace: need N*N)
1568: // *
1569: this._dgemm.Run("N", "N", M, N, N, ONE
1570: , A, offset_a, LDA, WORK, IU + o_work, LDWRKU, ZERO, ref U, offset_u
1571: , LDU);
1572: // *
1573: // * Copy right singular vectors of R to A
1574: // * (Workspace: need N*N)
1575: // *
1576: this._dlacpy.Run("F", N, N, WORK, IR + o_work, LDWRKR, ref A, offset_a
1577: , LDA);
1578: // *
1579: }
1580: else
1581: {
1582: // *
1583: // * Insufficient workspace for a fast algorithm
1584: // *
1585: ITAU = 1;
1586: IWORK = ITAU + N;
1587: // *
1588: // * Compute A=Q*R, copying result to U
1589: // * (Workspace: need 2*N, prefer N+N*NB)
1590: // *
1591: this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
1592: , LWORK - IWORK + 1, ref IERR);
1593: this._dlacpy.Run("L", M, N, A, offset_a, LDA, ref U, offset_u
1594: , LDU);
1595: // *
1596: // * Generate Q in U
1597: // * (Workspace: need 2*N, prefer N+N*NB)
1598: // *
1599: this._dorgqr.Run(M, N, N, ref U, offset_u, LDU, WORK, ITAU + o_work
1600: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1601: IE = ITAU;
1602: ITAUQ = IE + N;
1603: ITAUP = ITAUQ + N;
1604: IWORK = ITAUP + N;
1605: // *
1606: // * Zero out below R in A
1607: // *
1608: this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref A, 2+1 * LDA + o_a
1609: , LDA);
1610: // *
1611: // * Bidiagonalize R in A
1612: // * (Workspace: need 4*N, prefer 3*N+2*N*NB)
1613: // *
1614: this._dgebrd.Run(N, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
1615: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1616: // *
1617: // * Multiply Q in U by left vectors bidiagonalizing R
1618: // * (Workspace: need 3*N+M, prefer 3*N+M*NB)
1619: // *
1620: this._dormbr.Run("Q", "R", "N", M, N, N
1621: , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, IWORK + o_work
1622: , LWORK - IWORK + 1, ref IERR);
1623: // *
1624: // * Generate right vectors bidiagonalizing R in A
1625: // * (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
1626: // *
1627: this._dorgbr.Run("P", N, N, N, ref A, offset_a, LDA
1628: , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1629: IWORK = IE + N;
1630: // *
1631: // * Perform bidiagonal QR iteration, computing left
1632: // * singular vectors of A in U and computing right
1633: // * singular vectors of A in A
1634: // * (Workspace: need BDSPAC)
1635: // *
1636: this._dbdsqr.Run("U", N, N, M, 0, ref S, offset_s
1637: , ref WORK, IE + o_work, ref A, offset_a, LDA, ref U, offset_u, LDU, ref DUM, offset_dum
1638: , 1, ref WORK, IWORK + o_work, ref INFO);
1639: // *
1640: }
1641: // *
1642: }
1643: else
1644: {
1645: if (WNTVAS)
1646: {
1647: // *
1648: // * Path 6 (M much larger than N, JOBU='S', JOBVT='S'
1649: // * or 'A')
1650: // * N left singular vectors to be computed in U and
1651: // * N right singular vectors to be computed in VT
1652: // *
1653: if (LWORK >= N * N + Math.Max(4 * N, BDSPAC))
1654: {
1655: // *
1656: // * Sufficient workspace for a fast algorithm
1657: // *
1658: IU = 1;
1659: if (LWORK >= WRKBL + LDA * N)
1660: {
1661: // *
1662: // * WORK(IU) is LDA by N
1663: // *
1664: LDWRKU = LDA;
1665: }
1666: else
1667: {
1668: // *
1669: // * WORK(IU) is N by N
1670: // *
1671: LDWRKU = N;
1672: }
1673: ITAU = IU + LDWRKU * N;
1674: IWORK = ITAU + N;
1675: // *
1676: // * Compute A=Q*R
1677: // * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1678: // *
1679: this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
1680: , LWORK - IWORK + 1, ref IERR);
1681: // *
1682: // * Copy R to WORK(IU), zeroing out below it
1683: // *
1684: this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref WORK, IU + o_work
1685: , LDWRKU);
1686: this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref WORK, IU + 1 + o_work
1687: , LDWRKU);
1688: // *
1689: // * Generate Q in A
1690: // * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1691: // *
1692: this._dorgqr.Run(M, N, N, ref A, offset_a, LDA, WORK, ITAU + o_work
1693: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1694: IE = ITAU;
1695: ITAUQ = IE + N;
1696: ITAUP = ITAUQ + N;
1697: IWORK = ITAUP + N;
1698: // *
1699: // * Bidiagonalize R in WORK(IU), copying result to VT
1700: // * (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
1701: // *
1702: this._dgebrd.Run(N, N, ref WORK, IU + o_work, LDWRKU, ref S, offset_s, ref WORK, IE + o_work
1703: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1704: this._dlacpy.Run("U", N, N, WORK, IU + o_work, LDWRKU, ref VT, offset_vt
1705: , LDVT);
1706: // *
1707: // * Generate left bidiagonalizing vectors in WORK(IU)
1708: // * (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
1709: // *
1710: this._dorgbr.Run("Q", N, N, N, ref WORK, IU + o_work, LDWRKU
1711: , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1712: // *
1713: // * Generate right bidiagonalizing vectors in VT
1714: // * (Workspace: need N*N+4*N-1,
1715: // * prefer N*N+3*N+(N-1)*NB)
1716: // *
1717: this._dorgbr.Run("P", N, N, N, ref VT, offset_vt, LDVT
1718: , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1719: IWORK = IE + N;
1720: // *
1721: // * Perform bidiagonal QR iteration, computing left
1722: // * singular vectors of R in WORK(IU) and computing
1723: // * right singular vectors of R in VT
1724: // * (Workspace: need N*N+BDSPAC)
1725: // *
1726: this._dbdsqr.Run("U", N, N, N, 0, ref S, offset_s
1727: , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref WORK, IU + o_work, LDWRKU, ref DUM, offset_dum
1728: , 1, ref WORK, IWORK + o_work, ref INFO);
1729: // *
1730: // * Multiply Q in A by left singular vectors of R in
1731: // * WORK(IU), storing result in U
1732: // * (Workspace: need N*N)
1733: // *
1734: this._dgemm.Run("N", "N", M, N, N, ONE
1735: , A, offset_a, LDA, WORK, IU + o_work, LDWRKU, ZERO, ref U, offset_u
1736: , LDU);
1737: // *
1738: }
1739: else
1740: {
1741: // *
1742: // * Insufficient workspace for a fast algorithm
1743: // *
1744: ITAU = 1;
1745: IWORK = ITAU + N;
1746: // *
1747: // * Compute A=Q*R, copying result to U
1748: // * (Workspace: need 2*N, prefer N+N*NB)
1749: // *
1750: this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
1751: , LWORK - IWORK + 1, ref IERR);
1752: this._dlacpy.Run("L", M, N, A, offset_a, LDA, ref U, offset_u
1753: , LDU);
1754: // *
1755: // * Generate Q in U
1756: // * (Workspace: need 2*N, prefer N+N*NB)
1757: // *
1758: this._dorgqr.Run(M, N, N, ref U, offset_u, LDU, WORK, ITAU + o_work
1759: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1760: // *
1761: // * Copy R to VT, zeroing out below it
1762: // *
1763: this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref VT, offset_vt
1764: , LDVT);
1765: if (N > 1)
1766: {
1767: this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref VT, 2+1 * LDVT + o_vt
1768: , LDVT);
1769: }
1770: IE = ITAU;
1771: ITAUQ = IE + N;
1772: ITAUP = ITAUQ + N;
1773: IWORK = ITAUP + N;
1774: // *
1775: // * Bidiagonalize R in VT
1776: // * (Workspace: need 4*N, prefer 3*N+2*N*NB)
1777: // *
1778: this._dgebrd.Run(N, N, ref VT, offset_vt, LDVT, ref S, offset_s, ref WORK, IE + o_work
1779: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1780: // *
1781: // * Multiply Q in U by left bidiagonalizing vectors
1782: // * in VT
1783: // * (Workspace: need 3*N+M, prefer 3*N+M*NB)
1784: // *
1785: this._dormbr.Run("Q", "R", "N", M, N, N
1786: , ref VT, offset_vt, LDVT, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, IWORK + o_work
1787: , LWORK - IWORK + 1, ref IERR);
1788: // *
1789: // * Generate right bidiagonalizing vectors in VT
1790: // * (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
1791: // *
1792: this._dorgbr.Run("P", N, N, N, ref VT, offset_vt, LDVT
1793: , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1794: IWORK = IE + N;
1795: // *
1796: // * Perform bidiagonal QR iteration, computing left
1797: // * singular vectors of A in U and computing right
1798: // * singular vectors of A in VT
1799: // * (Workspace: need BDSPAC)
1800: // *
1801: this._dbdsqr.Run("U", N, N, M, 0, ref S, offset_s
1802: , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref U, offset_u, LDU, ref DUM, offset_dum
1803: , 1, ref WORK, IWORK + o_work, ref INFO);
1804: // *
1805: }
1806: // *
1807: }
1808: }
1809: }
1810: // *
1811: }
1812: else
1813: {
1814: if (WNTUA)
1815: {
1816: // *
1817: if (WNTVN)
1818: {
1819: // *
1820: // * Path 7 (M much larger than N, JOBU='A', JOBVT='N')
1821: // * M left singular vectors to be computed in U and
1822: // * no right singular vectors to be computed
1823: // *
1824: if (LWORK >= N * N + Math.Max(N + M, Math.Max(4 * N, BDSPAC)))
1825: {
1826: // *
1827: // * Sufficient workspace for a fast algorithm
1828: // *
1829: IR = 1;
1830: if (LWORK >= WRKBL + LDA * N)
1831: {
1832: // *
1833: // * WORK(IR) is LDA by N
1834: // *
1835: LDWRKR = LDA;
1836: }
1837: else
1838: {
1839: // *
1840: // * WORK(IR) is N by N
1841: // *
1842: LDWRKR = N;
1843: }
1844: ITAU = IR + LDWRKR * N;
1845: IWORK = ITAU + N;
1846: // *
1847: // * Compute A=Q*R, copying result to U
1848: // * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1849: // *
1850: this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
1851: , LWORK - IWORK + 1, ref IERR);
1852: this._dlacpy.Run("L", M, N, A, offset_a, LDA, ref U, offset_u
1853: , LDU);
1854: // *
1855: // * Copy R to WORK(IR), zeroing out below it
1856: // *
1857: this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref WORK, IR + o_work
1858: , LDWRKR);
1859: this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref WORK, IR + 1 + o_work
1860: , LDWRKR);
1861: // *
1862: // * Generate Q in U
1863: // * (Workspace: need N*N+N+M, prefer N*N+N+M*NB)
1864: // *
1865: this._dorgqr.Run(M, M, N, ref U, offset_u, LDU, WORK, ITAU + o_work
1866: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1867: IE = ITAU;
1868: ITAUQ = IE + N;
1869: ITAUP = ITAUQ + N;
1870: IWORK = ITAUP + N;
1871: // *
1872: // * Bidiagonalize R in WORK(IR)
1873: // * (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
1874: // *
1875: this._dgebrd.Run(N, N, ref WORK, IR + o_work, LDWRKR, ref S, offset_s, ref WORK, IE + o_work
1876: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1877: // *
1878: // * Generate left bidiagonalizing vectors in WORK(IR)
1879: // * (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
1880: // *
1881: this._dorgbr.Run("Q", N, N, N, ref WORK, IR + o_work, LDWRKR
1882: , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1883: IWORK = IE + N;
1884: // *
1885: // * Perform bidiagonal QR iteration, computing left
1886: // * singular vectors of R in WORK(IR)
1887: // * (Workspace: need N*N+BDSPAC)
1888: // *
1889: this._dbdsqr.Run("U", N, 0, N, 0, ref S, offset_s
1890: , ref WORK, IE + o_work, ref DUM, offset_dum, 1, ref WORK, IR + o_work, LDWRKR, ref DUM, offset_dum
1891: , 1, ref WORK, IWORK + o_work, ref INFO);
1892: // *
1893: // * Multiply Q in U by left singular vectors of R in
1894: // * WORK(IR), storing result in A
1895: // * (Workspace: need N*N)
1896: // *
1897: this._dgemm.Run("N", "N", M, N, N, ONE
1898: , U, offset_u, LDU, WORK, IR + o_work, LDWRKR, ZERO, ref A, offset_a
1899: , LDA);
1900: // *
1901: // * Copy left singular vectors of A from A to U
1902: // *
1903: this._dlacpy.Run("F", M, N, A, offset_a, LDA, ref U, offset_u
1904: , LDU);
1905: // *
1906: }
1907: else
1908: {
1909: // *
1910: // * Insufficient workspace for a fast algorithm
1911: // *
1912: ITAU = 1;
1913: IWORK = ITAU + N;
1914: // *
1915: // * Compute A=Q*R, copying result to U
1916: // * (Workspace: need 2*N, prefer N+N*NB)
1917: // *
1918: this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
1919: , LWORK - IWORK + 1, ref IERR);
1920: this._dlacpy.Run("L", M, N, A, offset_a, LDA, ref U, offset_u
1921: , LDU);
1922: // *
1923: // * Generate Q in U
1924: // * (Workspace: need N+M, prefer N+M*NB)
1925: // *
1926: this._dorgqr.Run(M, M, N, ref U, offset_u, LDU, WORK, ITAU + o_work
1927: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1928: IE = ITAU;
1929: ITAUQ = IE + N;
1930: ITAUP = ITAUQ + N;
1931: IWORK = ITAUP + N;
1932: // *
1933: // * Zero out below R in A
1934: // *
1935: this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref A, 2+1 * LDA + o_a
1936: , LDA);
1937: // *
1938: // * Bidiagonalize R in A
1939: // * (Workspace: need 4*N, prefer 3*N+2*N*NB)
1940: // *
1941: this._dgebrd.Run(N, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
1942: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1943: // *
1944: // * Multiply Q in U by left bidiagonalizing vectors
1945: // * in A
1946: // * (Workspace: need 3*N+M, prefer 3*N+M*NB)
1947: // *
1948: this._dormbr.Run("Q", "R", "N", M, N, N
1949: , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, IWORK + o_work
1950: , LWORK - IWORK + 1, ref IERR);
1951: IWORK = IE + N;
1952: // *
1953: // * Perform bidiagonal QR iteration, computing left
1954: // * singular vectors of A in U
1955: // * (Workspace: need BDSPAC)
1956: // *
1957: this._dbdsqr.Run("U", N, 0, M, 0, ref S, offset_s
1958: , ref WORK, IE + o_work, ref DUM, offset_dum, 1, ref U, offset_u, LDU, ref DUM, offset_dum
1959: , 1, ref WORK, IWORK + o_work, ref INFO);
1960: // *
1961: }
1962: // *
1963: }
1964: else
1965: {
1966: if (WNTVO)
1967: {
1968: // *
1969: // * Path 8 (M much larger than N, JOBU='A', JOBVT='O')
1970: // * M left singular vectors to be computed in U and
1971: // * N right singular vectors to be overwritten on A
1972: // *
1973: if (LWORK >= 2 * N * N + Math.Max(N + M, Math.Max(4 * N, BDSPAC)))
1974: {
1975: // *
1976: // * Sufficient workspace for a fast algorithm
1977: // *
1978: IU = 1;
1979: if (LWORK >= WRKBL + 2 * LDA * N)
1980: {
1981: // *
1982: // * WORK(IU) is LDA by N and WORK(IR) is LDA by N
1983: // *
1984: LDWRKU = LDA;
1985: IR = IU + LDWRKU * N;
1986: LDWRKR = LDA;
1987: }
1988: else
1989: {
1990: if (LWORK >= WRKBL + (LDA + N) * N)
1991: {
1992: // *
1993: // * WORK(IU) is LDA by N and WORK(IR) is N by N
1994: // *
1995: LDWRKU = LDA;
1996: IR = IU + LDWRKU * N;
1997: LDWRKR = N;
1998: }
1999: else
2000: {
2001: // *
2002: // * WORK(IU) is N by N and WORK(IR) is N by N
2003: // *
2004: LDWRKU = N;
2005: IR = IU + LDWRKU * N;
2006: LDWRKR = N;
2007: }
2008: }
2009: ITAU = IR + LDWRKR * N;
2010: IWORK = ITAU + N;
2011: // *
2012: // * Compute A=Q*R, copying result to U
2013: // * (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
2014: // *
2015: this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
2016: , LWORK - IWORK + 1, ref IERR);
2017: this._dlacpy.Run("L", M, N, A, offset_a, LDA, ref U, offset_u
2018: , LDU);
2019: // *
2020: // * Generate Q in U
2021: // * (Workspace: need 2*N*N+N+M, prefer 2*N*N+N+M*NB)
2022: // *
2023: this._dorgqr.Run(M, M, N, ref U, offset_u, LDU, WORK, ITAU + o_work
2024: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2025: // *
2026: // * Copy R to WORK(IU), zeroing out below it
2027: // *
2028: this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref WORK, IU + o_work
2029: , LDWRKU);
2030: this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref WORK, IU + 1 + o_work
2031: , LDWRKU);
2032: IE = ITAU;
2033: ITAUQ = IE + N;
2034: ITAUP = ITAUQ + N;
2035: IWORK = ITAUP + N;
2036: // *
2037: // * Bidiagonalize R in WORK(IU), copying result to
2038: // * WORK(IR)
2039: // * (Workspace: need 2*N*N+4*N,
2040: // * prefer 2*N*N+3*N+2*N*NB)
2041: // *
2042: this._dgebrd.Run(N, N, ref WORK, IU + o_work, LDWRKU, ref S, offset_s, ref WORK, IE + o_work
2043: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2044: this._dlacpy.Run("U", N, N, WORK, IU + o_work, LDWRKU, ref WORK, IR + o_work
2045: , LDWRKR);
2046: // *
2047: // * Generate left bidiagonalizing vectors in WORK(IU)
2048: // * (Workspace: need 2*N*N+4*N, prefer 2*N*N+3*N+N*NB)
2049: // *
2050: this._dorgbr.Run("Q", N, N, N, ref WORK, IU + o_work, LDWRKU
2051: , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2052: // *
2053: // * Generate right bidiagonalizing vectors in WORK(IR)
2054: // * (Workspace: need 2*N*N+4*N-1,
2055: // * prefer 2*N*N+3*N+(N-1)*NB)
2056: // *
2057: this._dorgbr.Run("P", N, N, N, ref WORK, IR + o_work, LDWRKR
2058: , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2059: IWORK = IE + N;
2060: // *
2061: // * Perform bidiagonal QR iteration, computing left
2062: // * singular vectors of R in WORK(IU) and computing
2063: // * right singular vectors of R in WORK(IR)
2064: // * (Workspace: need 2*N*N+BDSPAC)
2065: // *
2066: this._dbdsqr.Run("U", N, N, N, 0, ref S, offset_s
2067: , ref WORK, IE + o_work, ref WORK, IR + o_work, LDWRKR, ref WORK, IU + o_work, LDWRKU, ref DUM, offset_dum
2068: , 1, ref WORK, IWORK + o_work, ref INFO);
2069: // *
2070: // * Multiply Q in U by left singular vectors of R in
2071: // * WORK(IU), storing result in A
2072: // * (Workspace: need N*N)
2073: // *
2074: this._dgemm.Run("N", "N", M, N, N, ONE
2075: , U, offset_u, LDU, WORK, IU + o_work, LDWRKU, ZERO, ref A, offset_a
2076: , LDA);
2077: // *
2078: // * Copy left singular vectors of A from A to U
2079: // *
2080: this._dlacpy.Run("F", M, N, A, offset_a, LDA, ref U, offset_u
2081: , LDU);
2082: // *
2083: // * Copy right singular vectors of R from WORK(IR) to A
2084: // *
2085: this._dlacpy.Run("F", N, N, WORK, IR + o_work, LDWRKR, ref A, offset_a
2086: , LDA);
2087: // *
2088: }
2089: else
2090: {
2091: // *
2092: // * Insufficient workspace for a fast algorithm
2093: // *
2094: ITAU = 1;
2095: IWORK = ITAU + N;
2096: // *
2097: // * Compute A=Q*R, copying result to U
2098: // * (Workspace: need 2*N, prefer N+N*NB)
2099: // *
2100: this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
2101: , LWORK - IWORK + 1, ref IERR);
2102: this._dlacpy.Run("L", M, N, A, offset_a, LDA, ref U, offset_u
2103: , LDU);
2104: // *
2105: // * Generate Q in U
2106: // * (Workspace: need N+M, prefer N+M*NB)
2107: // *
2108: this._dorgqr.Run(M, M, N, ref U, offset_u, LDU, WORK, ITAU + o_work
2109: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2110: IE = ITAU;
2111: ITAUQ = IE + N;
2112: ITAUP = ITAUQ + N;
2113: IWORK = ITAUP + N;
2114: // *
2115: // * Zero out below R in A
2116: // *
2117: this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref A, 2+1 * LDA + o_a
2118: , LDA);
2119: // *
2120: // * Bidiagonalize R in A
2121: // * (Workspace: need 4*N, prefer 3*N+2*N*NB)
2122: // *
2123: this._dgebrd.Run(N, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
2124: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2125: // *
2126: // * Multiply Q in U by left bidiagonalizing vectors
2127: // * in A
2128: // * (Workspace: need 3*N+M, prefer 3*N+M*NB)
2129: // *
2130: this._dormbr.Run("Q", "R", "N", M, N, N
2131: , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, IWORK + o_work
2132: , LWORK - IWORK + 1, ref IERR);
2133: // *
2134: // * Generate right bidiagonalizing vectors in A
2135: // * (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
2136: // *
2137: this._dorgbr.Run("P", N, N, N, ref A, offset_a, LDA
2138: , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2139: IWORK = IE + N;
2140: // *
2141: // * Perform bidiagonal QR iteration, computing left
2142: // * singular vectors of A in U and computing right
2143: // * singular vectors of A in A
2144: // * (Workspace: need BDSPAC)
2145: // *
2146: this._dbdsqr.Run("U", N, N, M, 0, ref S, offset_s
2147: , ref WORK, IE + o_work, ref A, offset_a, LDA, ref U, offset_u, LDU, ref DUM, offset_dum
2148: , 1, ref WORK, IWORK + o_work, ref INFO);
2149: // *
2150: }
2151: // *
2152: }
2153: else
2154: {
2155: if (WNTVAS)
2156: {
2157: // *
2158: // * Path 9 (M much larger than N, JOBU='A', JOBVT='S'
2159: // * or 'A')
2160: // * M left singular vectors to be computed in U and
2161: // * N right singular vectors to be computed in VT
2162: // *
2163: if (LWORK >= N * N + Math.Max(N + M, Math.Max(4 * N, BDSPAC)))
2164: {
2165: // *
2166: // * Sufficient workspace for a fast algorithm
2167: // *
2168: IU = 1;
2169: if (LWORK >= WRKBL + LDA * N)
2170: {
2171: // *
2172: // * WORK(IU) is LDA by N
2173: // *
2174: LDWRKU = LDA;
2175: }
2176: else
2177: {
2178: // *
2179: // * WORK(IU) is N by N
2180: // *
2181: LDWRKU = N;
2182: }
2183: ITAU = IU + LDWRKU * N;
2184: IWORK = ITAU + N;
2185: // *
2186: // * Compute A=Q*R, copying result to U
2187: // * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
2188: // *
2189: this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
2190: , LWORK - IWORK + 1, ref IERR);
2191: this._dlacpy.Run("L", M, N, A, offset_a, LDA, ref U, offset_u
2192: , LDU);
2193: // *
2194: // * Generate Q in U
2195: // * (Workspace: need N*N+N+M, prefer N*N+N+M*NB)
2196: // *
2197: this._dorgqr.Run(M, M, N, ref U, offset_u, LDU, WORK, ITAU + o_work
2198: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2199: // *
2200: // * Copy R to WORK(IU), zeroing out below it
2201: // *
2202: this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref WORK, IU + o_work
2203: , LDWRKU);
2204: this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref WORK, IU + 1 + o_work
2205: , LDWRKU);
2206: IE = ITAU;
2207: ITAUQ = IE + N;
2208: ITAUP = ITAUQ + N;
2209: IWORK = ITAUP + N;
2210: // *
2211: // * Bidiagonalize R in WORK(IU), copying result to VT
2212: // * (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
2213: // *
2214: this._dgebrd.Run(N, N, ref WORK, IU + o_work, LDWRKU, ref S, offset_s, ref WORK, IE + o_work
2215: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2216: this._dlacpy.Run("U", N, N, WORK, IU + o_work, LDWRKU, ref VT, offset_vt
2217: , LDVT);
2218: // *
2219: // * Generate left bidiagonalizing vectors in WORK(IU)
2220: // * (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
2221: // *
2222: this._dorgbr.Run("Q", N, N, N, ref WORK, IU + o_work, LDWRKU
2223: , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2224: // *
2225: // * Generate right bidiagonalizing vectors in VT
2226: // * (Workspace: need N*N+4*N-1,
2227: // * prefer N*N+3*N+(N-1)*NB)
2228: // *
2229: this._dorgbr.Run("P", N, N, N, ref VT, offset_vt, LDVT
2230: , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2231: IWORK = IE + N;
2232: // *
2233: // * Perform bidiagonal QR iteration, computing left
2234: // * singular vectors of R in WORK(IU) and computing
2235: // * right singular vectors of R in VT
2236: // * (Workspace: need N*N+BDSPAC)
2237: // *
2238: this._dbdsqr.Run("U", N, N, N, 0, ref S, offset_s
2239: , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref WORK, IU + o_work, LDWRKU, ref DUM, offset_dum
2240: , 1, ref WORK, IWORK + o_work, ref INFO);
2241: // *
2242: // * Multiply Q in U by left singular vectors of R in
2243: // * WORK(IU), storing result in A
2244: // * (Workspace: need N*N)
2245: // *
2246: this._dgemm.Run("N", "N", M, N, N, ONE
2247: , U, offset_u, LDU, WORK, IU + o_work, LDWRKU, ZERO, ref A, offset_a
2248: , LDA);
2249: // *
2250: // * Copy left singular vectors of A from A to U
2251: // *
2252: this._dlacpy.Run("F", M, N, A, offset_a, LDA, ref U, offset_u
2253: , LDU);
2254: // *
2255: }
2256: else
2257: {
2258: // *
2259: // * Insufficient workspace for a fast algorithm
2260: // *
2261: ITAU = 1;
2262: IWORK = ITAU + N;
2263: // *
2264: // * Compute A=Q*R, copying result to U
2265: // * (Workspace: need 2*N, prefer N+N*NB)
2266: // *
2267: this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
2268: , LWORK - IWORK + 1, ref IERR);
2269: this._dlacpy.Run("L", M, N, A, offset_a, LDA, ref U, offset_u
2270: , LDU);
2271: // *
2272: // * Generate Q in U
2273: // * (Workspace: need N+M, prefer N+M*NB)
2274: // *
2275: this._dorgqr.Run(M, M, N, ref U, offset_u, LDU, WORK, ITAU + o_work
2276: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2277: // *
2278: // * Copy R from A to VT, zeroing out below it
2279: // *
2280: this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref VT, offset_vt
2281: , LDVT);
2282: if (N > 1)
2283: {
2284: this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref VT, 2+1 * LDVT + o_vt
2285: , LDVT);
2286: }
2287: IE = ITAU;
2288: ITAUQ = IE + N;
2289: ITAUP = ITAUQ + N;
2290: IWORK = ITAUP + N;
2291: // *
2292: // * Bidiagonalize R in VT
2293: // * (Workspace: need 4*N, prefer 3*N+2*N*NB)
2294: // *
2295: this._dgebrd.Run(N, N, ref VT, offset_vt, LDVT, ref S, offset_s, ref WORK, IE + o_work
2296: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2297: // *
2298: // * Multiply Q in U by left bidiagonalizing vectors
2299: // * in VT
2300: // * (Workspace: need 3*N+M, prefer 3*N+M*NB)
2301: // *
2302: this._dormbr.Run("Q", "R", "N", M, N, N
2303: , ref VT, offset_vt, LDVT, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, IWORK + o_work
2304: , LWORK - IWORK + 1, ref IERR);
2305: // *
2306: // * Generate right bidiagonalizing vectors in VT
2307: // * (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
2308: // *
2309: this._dorgbr.Run("P", N, N, N, ref VT, offset_vt, LDVT
2310: , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2311: IWORK = IE + N;
2312: // *
2313: // * Perform bidiagonal QR iteration, computing left
2314: // * singular vectors of A in U and computing right
2315: // * singular vectors of A in VT
2316: // * (Workspace: need BDSPAC)
2317: // *
2318: this._dbdsqr.Run("U", N, N, M, 0, ref S, offset_s
2319: , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref U, offset_u, LDU, ref DUM, offset_dum
2320: , 1, ref WORK, IWORK + o_work, ref INFO);
2321: // *
2322: }
2323: // *
2324: }
2325: }
2326: }
2327: // *
2328: }
2329: }
2330: }
2331: }
2332: }
2333: // *
2334: }
2335: else
2336: {
2337: // *
2338: // * M .LT. MNTHR
2339: // *
2340: // * Path 10 (M at least N, but not much larger)
2341: // * Reduce to bidiagonal form without QR decomposition
2342: // *
2343: IE = 1;
2344: ITAUQ = IE + N;
2345: ITAUP = ITAUQ + N;
2346: IWORK = ITAUP + N;
2347: // *
2348: // * Bidiagonalize A
2349: // * (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB)
2350: // *
2351: this._dgebrd.Run(M, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
2352: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2353: if (WNTUAS)
2354: {
2355: // *
2356: // * If left singular vectors desired in U, copy result to U
2357: // * and generate left bidiagonalizing vectors in U
2358: // * (Workspace: need 3*N+NCU, prefer 3*N+NCU*NB)
2359: // *
2360: this._dlacpy.Run("L", M, N, A, offset_a, LDA, ref U, offset_u
2361: , LDU);
2362: if (WNTUS) NCU = N;
2363: if (WNTUA) NCU = M;
2364: this._dorgbr.Run("Q", M, NCU, N, ref U, offset_u, LDU
2365: , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2366: }
2367: if (WNTVAS)
2368: {
2369: // *
2370: // * If right singular vectors desired in VT, copy result to
2371: // * VT and generate right bidiagonalizing vectors in VT
2372: // * (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
2373: // *
2374: this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref VT, offset_vt
2375: , LDVT);
2376: this._dorgbr.Run("P", N, N, N, ref VT, offset_vt, LDVT
2377: , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2378: }
2379: if (WNTUO)
2380: {
2381: // *
2382: // * If left singular vectors desired in A, generate left
2383: // * bidiagonalizing vectors in A
2384: // * (Workspace: need 4*N, prefer 3*N+N*NB)
2385: // *
2386: this._dorgbr.Run("Q", M, N, N, ref A, offset_a, LDA
2387: , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2388: }
2389: if (WNTVO)
2390: {
2391: // *
2392: // * If right singular vectors desired in A, generate right
2393: // * bidiagonalizing vectors in A
2394: // * (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
2395: // *
2396: this._dorgbr.Run("P", N, N, N, ref A, offset_a, LDA
2397: , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2398: }
2399: IWORK = IE + N;
2400: if (WNTUAS || WNTUO) NRU = M;
2401: if (WNTUN) NRU = 0;
2402: if (WNTVAS || WNTVO) NCVT = N;
2403: if (WNTVN) NCVT = 0;
2404: if ((!WNTUO) && (!WNTVO))
2405: {
2406: // *
2407: // * Perform bidiagonal QR iteration, if desired, computing
2408: // * left singular vectors in U and computing right singular
2409: // * vectors in VT
2410: // * (Workspace: need BDSPAC)
2411: // *
2412: this._dbdsqr.Run("U", N, NCVT, NRU, 0, ref S, offset_s
2413: , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref U, offset_u, LDU, ref DUM, offset_dum
2414: , 1, ref WORK, IWORK + o_work, ref INFO);
2415: }
2416: else
2417: {
2418: if ((!WNTUO) && WNTVO)
2419: {
2420: // *
2421: // * Perform bidiagonal QR iteration, if desired, computing
2422: // * left singular vectors in U and computing right singular
2423: // * vectors in A
2424: // * (Workspace: need BDSPAC)
2425: // *
2426: this._dbdsqr.Run("U", N, NCVT, NRU, 0, ref S, offset_s
2427: , ref WORK, IE + o_work, ref A, offset_a, LDA, ref U, offset_u, LDU, ref DUM, offset_dum
2428: , 1, ref WORK, IWORK + o_work, ref INFO);
2429: }
2430: else
2431: {
2432: // *
2433: // * Perform bidiagonal QR iteration, if desired, computing
2434: // * left singular vectors in A and computing right singular
2435: // * vectors in VT
2436: // * (Workspace: need BDSPAC)
2437: // *
2438: this._dbdsqr.Run("U", N, NCVT, NRU, 0, ref S, offset_s
2439: , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref A, offset_a, LDA, ref DUM, offset_dum
2440: , 1, ref WORK, IWORK + o_work, ref INFO);
2441: }
2442: }
2443: // *
2444: }
2445: // *
2446: }
2447: else
2448: {
2449: // *
2450: // * A has more columns than rows. If A has sufficiently more
2451: // * columns than rows, first reduce using the LQ decomposition (if
2452: // * sufficient workspace available)
2453: // *
2454: if (N >= MNTHR)
2455: {
2456: // *
2457: if (WNTVN)
2458: {
2459: // *
2460: // * Path 1t(N much larger than M, JOBVT='N')
2461: // * No right singular vectors to be computed
2462: // *
2463: ITAU = 1;
2464: IWORK = ITAU + M;
2465: // *
2466: // * Compute A=L*Q
2467: // * (Workspace: need 2*M, prefer M+M*NB)
2468: // *
2469: this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
2470: , LWORK - IWORK + 1, ref IERR);
2471: // *
2472: // * Zero out above L
2473: // *
2474: this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref A, 1+2 * LDA + o_a
2475: , LDA);
2476: IE = 1;
2477: ITAUQ = IE + M;
2478: ITAUP = ITAUQ + M;
2479: IWORK = ITAUP + M;
2480: // *
2481: // * Bidiagonalize L in A
2482: // * (Workspace: need 4*M, prefer 3*M+2*M*NB)
2483: // *
2484: this._dgebrd.Run(M, M, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
2485: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2486: if (WNTUO || WNTUAS)
2487: {
2488: // *
2489: // * If left singular vectors desired, generate Q
2490: // * (Workspace: need 4*M, prefer 3*M+M*NB)
2491: // *
2492: this._dorgbr.Run("Q", M, M, M, ref A, offset_a, LDA
2493: , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2494: }
2495: IWORK = IE + M;
2496: NRU = 0;
2497: if (WNTUO || WNTUAS) NRU = M;
2498: // *
2499: // * Perform bidiagonal QR iteration, computing left singular
2500: // * vectors of A in A if desired
2501: // * (Workspace: need BDSPAC)
2502: // *
2503: this._dbdsqr.Run("U", M, 0, NRU, 0, ref S, offset_s
2504: , ref WORK, IE + o_work, ref DUM, offset_dum, 1, ref A, offset_a, LDA, ref DUM, offset_dum
2505: , 1, ref WORK, IWORK + o_work, ref INFO);
2506: // *
2507: // * If left singular vectors desired in U, copy them there
2508: // *
2509: if (WNTUAS)
2510: {
2511: this._dlacpy.Run("F", M, M, A, offset_a, LDA, ref U, offset_u
2512: , LDU);
2513: }
2514: // *
2515: }
2516: else
2517: {
2518: if (WNTVO && WNTUN)
2519: {
2520: // *
2521: // * Path 2t(N much larger than M, JOBU='N', JOBVT='O')
2522: // * M right singular vectors to be overwritten on A and
2523: // * no left singular vectors to be computed
2524: // *
2525: if (LWORK >= M * M + Math.Max(4 * M, BDSPAC))
2526: {
2527: // *
2528: // * Sufficient workspace for a fast algorithm
2529: // *
2530: IR = 1;
2531: if (LWORK >= Math.Max(WRKBL, LDA * N + M) + LDA * M)
2532: {
2533: // *
2534: // * WORK(IU) is LDA by N and WORK(IR) is LDA by M
2535: // *
2536: LDWRKU = LDA;
2537: CHUNK = N;
2538: LDWRKR = LDA;
2539: }
2540: else
2541: {
2542: if (LWORK >= Math.Max(WRKBL, LDA * N + M) + M * M)
2543: {
2544: // *
2545: // * WORK(IU) is LDA by N and WORK(IR) is M by M
2546: // *
2547: LDWRKU = LDA;
2548: CHUNK = N;
2549: LDWRKR = M;
2550: }
2551: else
2552: {
2553: // *
2554: // * WORK(IU) is M by CHUNK and WORK(IR) is M by M
2555: // *
2556: LDWRKU = M;
2557: CHUNK = (LWORK - M * M - M) / M;
2558: LDWRKR = M;
2559: }
2560: }
2561: ITAU = IR + LDWRKR * M;
2562: IWORK = ITAU + M;
2563: // *
2564: // * Compute A=L*Q
2565: // * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
2566: // *
2567: this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
2568: , LWORK - IWORK + 1, ref IERR);
2569: // *
2570: // * Copy L to WORK(IR) and zero out above it
2571: // *
2572: this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref WORK, IR + o_work
2573: , LDWRKR);
2574: this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref WORK, IR + LDWRKR + o_work
2575: , LDWRKR);
2576: // *
2577: // * Generate Q in A
2578: // * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
2579: // *
2580: this._dorglq.Run(M, N, M, ref A, offset_a, LDA, WORK, ITAU + o_work
2581: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2582: IE = ITAU;
2583: ITAUQ = IE + M;
2584: ITAUP = ITAUQ + M;
2585: IWORK = ITAUP + M;
2586: // *
2587: // * Bidiagonalize L in WORK(IR)
2588: // * (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
2589: // *
2590: this._dgebrd.Run(M, M, ref WORK, IR + o_work, LDWRKR, ref S, offset_s, ref WORK, IE + o_work
2591: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2592: // *
2593: // * Generate right vectors bidiagonalizing L
2594: // * (Workspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)
2595: // *
2596: this._dorgbr.Run("P", M, M, M, ref WORK, IR + o_work, LDWRKR
2597: , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2598: IWORK = IE + M;
2599: // *
2600: // * Perform bidiagonal QR iteration, computing right
2601: // * singular vectors of L in WORK(IR)
2602: // * (Workspace: need M*M+BDSPAC)
2603: // *
2604: this._dbdsqr.Run("U", M, M, 0, 0, ref S, offset_s
2605: , ref WORK, IE + o_work, ref WORK, IR + o_work, LDWRKR, ref DUM, offset_dum, 1, ref DUM, offset_dum
2606: , 1, ref WORK, IWORK + o_work, ref INFO);
2607: IU = IE + M;
2608: // *
2609: // * Multiply right singular vectors of L in WORK(IR) by Q
2610: // * in A, storing result in WORK(IU) and copying to A
2611: // * (Workspace: need M*M+2*M, prefer M*M+M*N+M)
2612: // *
2613: for (I = 1; (CHUNK >= 0) ? (I <= N) : (I >= N); I += CHUNK)
2614: {
2615: BLK = Math.Min(N - I + 1, CHUNK);
2616: this._dgemm.Run("N", "N", M, BLK, M, ONE
2617: , WORK, IR + o_work, LDWRKR, A, 1+I * LDA + o_a, LDA, ZERO, ref WORK, IU + o_work
2618: , LDWRKU);
2619: this._dlacpy.Run("F", M, BLK, WORK, IU + o_work, LDWRKU, ref A, 1+I * LDA + o_a
2620: , LDA);
2621: }
2622: // *
2623: }
2624: else
2625: {
2626: // *
2627: // * Insufficient workspace for a fast algorithm
2628: // *
2629: IE = 1;
2630: ITAUQ = IE + M;
2631: ITAUP = ITAUQ + M;
2632: IWORK = ITAUP + M;
2633: // *
2634: // * Bidiagonalize A
2635: // * (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
2636: // *
2637: this._dgebrd.Run(M, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
2638: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2639: // *
2640: // * Generate right vectors bidiagonalizing A
2641: // * (Workspace: need 4*M, prefer 3*M+M*NB)
2642: // *
2643: this._dorgbr.Run("P", M, N, M, ref A, offset_a, LDA
2644: , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2645: IWORK = IE + M;
2646: // *
2647: // * Perform bidiagonal QR iteration, computing right
2648: // * singular vectors of A in A
2649: // * (Workspace: need BDSPAC)
2650: // *
2651: this._dbdsqr.Run("L", M, N, 0, 0, ref S, offset_s
2652: , ref WORK, IE + o_work, ref A, offset_a, LDA, ref DUM, offset_dum, 1, ref DUM, offset_dum
2653: , 1, ref WORK, IWORK + o_work, ref INFO);
2654: // *
2655: }
2656: // *
2657: }
2658: else
2659: {
2660: if (WNTVO && WNTUAS)
2661: {
2662: // *
2663: // * Path 3t(N much larger than M, JOBU='S' or 'A', JOBVT='O')
2664: // * M right singular vectors to be overwritten on A and
2665: // * M left singular vectors to be computed in U
2666: // *
2667: if (LWORK >= M * M + Math.Max(4 * M, BDSPAC))
2668: {
2669: // *
2670: // * Sufficient workspace for a fast algorithm
2671: // *
2672: IR = 1;
2673: if (LWORK >= Math.Max(WRKBL, LDA * N + M) + LDA * M)
2674: {
2675: // *
2676: // * WORK(IU) is LDA by N and WORK(IR) is LDA by M
2677: // *
2678: LDWRKU = LDA;
2679: CHUNK = N;
2680: LDWRKR = LDA;
2681: }
2682: else
2683: {
2684: if (LWORK >= Math.Max(WRKBL, LDA * N + M) + M * M)
2685: {
2686: // *
2687: // * WORK(IU) is LDA by N and WORK(IR) is M by M
2688: // *
2689: LDWRKU = LDA;
2690: CHUNK = N;
2691: LDWRKR = M;
2692: }
2693: else
2694: {
2695: // *
2696: // * WORK(IU) is M by CHUNK and WORK(IR) is M by M
2697: // *
2698: LDWRKU = M;
2699: CHUNK = (LWORK - M * M - M) / M;
2700: LDWRKR = M;
2701: }
2702: }
2703: ITAU = IR + LDWRKR * M;
2704: IWORK = ITAU + M;
2705: // *
2706: // * Compute A=L*Q
2707: // * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
2708: // *
2709: this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
2710: , LWORK - IWORK + 1, ref IERR);
2711: // *
2712: // * Copy L to U, zeroing about above it
2713: // *
2714: this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref U, offset_u
2715: , LDU);
2716: this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref U, 1+2 * LDU + o_u
2717: , LDU);
2718: // *
2719: // * Generate Q in A
2720: // * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
2721: // *
2722: this._dorglq.Run(M, N, M, ref A, offset_a, LDA, WORK, ITAU + o_work
2723: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2724: IE = ITAU;
2725: ITAUQ = IE + M;
2726: ITAUP = ITAUQ + M;
2727: IWORK = ITAUP + M;
2728: // *
2729: // * Bidiagonalize L in U, copying result to WORK(IR)
2730: // * (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
2731: // *
2732: this._dgebrd.Run(M, M, ref U, offset_u, LDU, ref S, offset_s, ref WORK, IE + o_work
2733: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2734: this._dlacpy.Run("U", M, M, U, offset_u, LDU, ref WORK, IR + o_work
2735: , LDWRKR);
2736: // *
2737: // * Generate right vectors bidiagonalizing L in WORK(IR)
2738: // * (Workspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)
2739: // *
2740: this._dorgbr.Run("P", M, M, M, ref WORK, IR + o_work, LDWRKR
2741: , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2742: // *
2743: // * Generate left vectors bidiagonalizing L in U
2744: // * (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB)
2745: // *
2746: this._dorgbr.Run("Q", M, M, M, ref U, offset_u, LDU
2747: , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2748: IWORK = IE + M;
2749: // *
2750: // * Perform bidiagonal QR iteration, computing left
2751: // * singular vectors of L in U, and computing right
2752: // * singular vectors of L in WORK(IR)
2753: // * (Workspace: need M*M+BDSPAC)
2754: // *
2755: this._dbdsqr.Run("U", M, M, M, 0, ref S, offset_s
2756: , ref WORK, IE + o_work, ref WORK, IR + o_work, LDWRKR, ref U, offset_u, LDU, ref DUM, offset_dum
2757: , 1, ref WORK, IWORK + o_work, ref INFO);
2758: IU = IE + M;
2759: // *
2760: // * Multiply right singular vectors of L in WORK(IR) by Q
2761: // * in A, storing result in WORK(IU) and copying to A
2762: // * (Workspace: need M*M+2*M, prefer M*M+M*N+M))
2763: // *
2764: for (I = 1; (CHUNK >= 0) ? (I <= N) : (I >= N); I += CHUNK)
2765: {
2766: BLK = Math.Min(N - I + 1, CHUNK);
2767: this._dgemm.Run("N", "N", M, BLK, M, ONE
2768: , WORK, IR + o_work, LDWRKR, A, 1+I * LDA + o_a, LDA, ZERO, ref WORK, IU + o_work
2769: , LDWRKU);
2770: this._dlacpy.Run("F", M, BLK, WORK, IU + o_work, LDWRKU, ref A, 1+I * LDA + o_a
2771: , LDA);
2772: }
2773: // *
2774: }
2775: else
2776: {
2777: // *
2778: // * Insufficient workspace for a fast algorithm
2779: // *
2780: ITAU = 1;
2781: IWORK = ITAU + M;
2782: // *
2783: // * Compute A=L*Q
2784: // * (Workspace: need 2*M, prefer M+M*NB)
2785: // *
2786: this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
2787: , LWORK - IWORK + 1, ref IERR);
2788: // *
2789: // * Copy L to U, zeroing out above it
2790: // *
2791: this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref U, offset_u
2792: , LDU);
2793: this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref U, 1+2 * LDU + o_u
2794: , LDU);
2795: // *
2796: // * Generate Q in A
2797: // * (Workspace: need 2*M, prefer M+M*NB)
2798: // *
2799: this._dorglq.Run(M, N, M, ref A, offset_a, LDA, WORK, ITAU + o_work
2800: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2801: IE = ITAU;
2802: ITAUQ = IE + M;
2803: ITAUP = ITAUQ + M;
2804: IWORK = ITAUP + M;
2805: // *
2806: // * Bidiagonalize L in U
2807: // * (Workspace: need 4*M, prefer 3*M+2*M*NB)
2808: // *
2809: this._dgebrd.Run(M, M, ref U, offset_u, LDU, ref S, offset_s, ref WORK, IE + o_work
2810: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2811: // *
2812: // * Multiply right vectors bidiagonalizing L by Q in A
2813: // * (Workspace: need 3*M+N, prefer 3*M+N*NB)
2814: // *
2815: this._dormbr.Run("P", "L", "T", M, N, M
2816: , ref U, offset_u, LDU, WORK, ITAUP + o_work, ref A, offset_a, LDA, ref WORK, IWORK + o_work
2817: , LWORK - IWORK + 1, ref IERR);
2818: // *
2819: // * Generate left vectors bidiagonalizing L in U
2820: // * (Workspace: need 4*M, prefer 3*M+M*NB)
2821: // *
2822: this._dorgbr.Run("Q", M, M, M, ref U, offset_u, LDU
2823: , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2824: IWORK = IE + M;
2825: // *
2826: // * Perform bidiagonal QR iteration, computing left
2827: // * singular vectors of A in U and computing right
2828: // * singular vectors of A in A
2829: // * (Workspace: need BDSPAC)
2830: // *
2831: this._dbdsqr.Run("U", M, N, M, 0, ref S, offset_s
2832: , ref WORK, IE + o_work, ref A, offset_a, LDA, ref U, offset_u, LDU, ref DUM, offset_dum
2833: , 1, ref WORK, IWORK + o_work, ref INFO);
2834: // *
2835: }
2836: // *
2837: }
2838: else
2839: {
2840: if (WNTVS)
2841: {
2842: // *
2843: if (WNTUN)
2844: {
2845: // *
2846: // * Path 4t(N much larger than M, JOBU='N', JOBVT='S')
2847: // * M right singular vectors to be computed in VT and
2848: // * no left singular vectors to be computed
2849: // *
2850: if (LWORK >= M * M + Math.Max(4 * M, BDSPAC))
2851: {
2852: // *
2853: // * Sufficient workspace for a fast algorithm
2854: // *
2855: IR = 1;
2856: if (LWORK >= WRKBL + LDA * M)
2857: {
2858: // *
2859: // * WORK(IR) is LDA by M
2860: // *
2861: LDWRKR = LDA;
2862: }
2863: else
2864: {
2865: // *
2866: // * WORK(IR) is M by M
2867: // *
2868: LDWRKR = M;
2869: }
2870: ITAU = IR + LDWRKR * M;
2871: IWORK = ITAU + M;
2872: // *
2873: // * Compute A=L*Q
2874: // * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
2875: // *
2876: this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
2877: , LWORK - IWORK + 1, ref IERR);
2878: // *
2879: // * Copy L to WORK(IR), zeroing out above it
2880: // *
2881: this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref WORK, IR + o_work
2882: , LDWRKR);
2883: this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref WORK, IR + LDWRKR + o_work
2884: , LDWRKR);
2885: // *
2886: // * Generate Q in A
2887: // * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
2888: // *
2889: this._dorglq.Run(M, N, M, ref A, offset_a, LDA, WORK, ITAU + o_work
2890: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2891: IE = ITAU;
2892: ITAUQ = IE + M;
2893: ITAUP = ITAUQ + M;
2894: IWORK = ITAUP + M;
2895: // *
2896: // * Bidiagonalize L in WORK(IR)
2897: // * (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
2898: // *
2899: this._dgebrd.Run(M, M, ref WORK, IR + o_work, LDWRKR, ref S, offset_s, ref WORK, IE + o_work
2900: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2901: // *
2902: // * Generate right vectors bidiagonalizing L in
2903: // * WORK(IR)
2904: // * (Workspace: need M*M+4*M, prefer M*M+3*M+(M-1)*NB)
2905: // *
2906: this._dorgbr.Run("P", M, M, M, ref WORK, IR + o_work, LDWRKR
2907: , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2908: IWORK = IE + M;
2909: // *
2910: // * Perform bidiagonal QR iteration, computing right
2911: // * singular vectors of L in WORK(IR)
2912: // * (Workspace: need M*M+BDSPAC)
2913: // *
2914: this._dbdsqr.Run("U", M, M, 0, 0, ref S, offset_s
2915: , ref WORK, IE + o_work, ref WORK, IR + o_work, LDWRKR, ref DUM, offset_dum, 1, ref DUM, offset_dum
2916: , 1, ref WORK, IWORK + o_work, ref INFO);
2917: // *
2918: // * Multiply right singular vectors of L in WORK(IR) by
2919: // * Q in A, storing result in VT
2920: // * (Workspace: need M*M)
2921: // *
2922: this._dgemm.Run("N", "N", M, N, M, ONE
2923: , WORK, IR + o_work, LDWRKR, A, offset_a, LDA, ZERO, ref VT, offset_vt
2924: , LDVT);
2925: // *
2926: }
2927: else
2928: {
2929: // *
2930: // * Insufficient workspace for a fast algorithm
2931: // *
2932: ITAU = 1;
2933: IWORK = ITAU + M;
2934: // *
2935: // * Compute A=L*Q
2936: // * (Workspace: need 2*M, prefer M+M*NB)
2937: // *
2938: this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
2939: , LWORK - IWORK + 1, ref IERR);
2940: // *
2941: // * Copy result to VT
2942: // *
2943: this._dlacpy.Run("U", M, N, A, offset_a, LDA, ref VT, offset_vt
2944: , LDVT);
2945: // *
2946: // * Generate Q in VT
2947: // * (Workspace: need 2*M, prefer M+M*NB)
2948: // *
2949: this._dorglq.Run(M, N, M, ref VT, offset_vt, LDVT, WORK, ITAU + o_work
2950: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2951: IE = ITAU;
2952: ITAUQ = IE + M;
2953: ITAUP = ITAUQ + M;
2954: IWORK = ITAUP + M;
2955: // *
2956: // * Zero out above L in A
2957: // *
2958: this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref A, 1+2 * LDA + o_a
2959: , LDA);
2960: // *
2961: // * Bidiagonalize L in A
2962: // * (Workspace: need 4*M, prefer 3*M+2*M*NB)
2963: // *
2964: this._dgebrd.Run(M, M, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
2965: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2966: // *
2967: // * Multiply right vectors bidiagonalizing L by Q in VT
2968: // * (Workspace: need 3*M+N, prefer 3*M+N*NB)
2969: // *
2970: this._dormbr.Run("P", "L", "T", M, N, M
2971: , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, IWORK + o_work
2972: , LWORK - IWORK + 1, ref IERR);
2973: IWORK = IE + M;
2974: // *
2975: // * Perform bidiagonal QR iteration, computing right
2976: // * singular vectors of A in VT
2977: // * (Workspace: need BDSPAC)
2978: // *
2979: this._dbdsqr.Run("U", M, N, 0, 0, ref S, offset_s
2980: , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref DUM, offset_dum, 1, ref DUM, offset_dum
2981: , 1, ref WORK, IWORK + o_work, ref INFO);
2982: // *
2983: }
2984: // *
2985: }
2986: else
2987: {
2988: if (WNTUO)
2989: {
2990: // *
2991: // * Path 5t(N much larger than M, JOBU='O', JOBVT='S')
2992: // * M right singular vectors to be computed in VT and
2993: // * M left singular vectors to be overwritten on A
2994: // *
2995: if (LWORK >= 2 * M * M + Math.Max(4 * M, BDSPAC))
2996: {
2997: // *
2998: // * Sufficient workspace for a fast algorithm
2999: // *
3000: IU = 1;
3001: if (LWORK >= WRKBL + 2 * LDA * M)
3002: {
3003: // *
3004: // * WORK(IU) is LDA by M and WORK(IR) is LDA by M
3005: // *
3006: LDWRKU = LDA;
3007: IR = IU + LDWRKU * M;
3008: LDWRKR = LDA;
3009: }
3010: else
3011: {
3012: if (LWORK >= WRKBL + (LDA + M) * M)
3013: {
3014: // *
3015: // * WORK(IU) is LDA by M and WORK(IR) is M by M
3016: // *
3017: LDWRKU = LDA;
3018: IR = IU + LDWRKU * M;
3019: LDWRKR = M;
3020: }
3021: else
3022: {
3023: // *
3024: // * WORK(IU) is M by M and WORK(IR) is M by M
3025: // *
3026: LDWRKU = M;
3027: IR = IU + LDWRKU * M;
3028: LDWRKR = M;
3029: }
3030: }
3031: ITAU = IR + LDWRKR * M;
3032: IWORK = ITAU + M;
3033: // *
3034: // * Compute A=L*Q
3035: // * (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
3036: // *
3037: this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
3038: , LWORK - IWORK + 1, ref IERR);
3039: // *
3040: // * Copy L to WORK(IU), zeroing out below it
3041: // *
3042: this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref WORK, IU + o_work
3043: , LDWRKU);
3044: this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref WORK, IU + LDWRKU + o_work
3045: , LDWRKU);
3046: // *
3047: // * Generate Q in A
3048: // * (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
3049: // *
3050: this._dorglq.Run(M, N, M, ref A, offset_a, LDA, WORK, ITAU + o_work
3051: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3052: IE = ITAU;
3053: ITAUQ = IE + M;
3054: ITAUP = ITAUQ + M;
3055: IWORK = ITAUP + M;
3056: // *
3057: // * Bidiagonalize L in WORK(IU), copying result to
3058: // * WORK(IR)
3059: // * (Workspace: need 2*M*M+4*M,
3060: // * prefer 2*M*M+3*M+2*M*NB)
3061: // *
3062: this._dgebrd.Run(M, M, ref WORK, IU + o_work, LDWRKU, ref S, offset_s, ref WORK, IE + o_work
3063: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3064: this._dlacpy.Run("L", M, M, WORK, IU + o_work, LDWRKU, ref WORK, IR + o_work
3065: , LDWRKR);
3066: // *
3067: // * Generate right bidiagonalizing vectors in WORK(IU)
3068: // * (Workspace: need 2*M*M+4*M-1,
3069: // * prefer 2*M*M+3*M+(M-1)*NB)
3070: // *
3071: this._dorgbr.Run("P", M, M, M, ref WORK, IU + o_work, LDWRKU
3072: , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3073: // *
3074: // * Generate left bidiagonalizing vectors in WORK(IR)
3075: // * (Workspace: need 2*M*M+4*M, prefer 2*M*M+3*M+M*NB)
3076: // *
3077: this._dorgbr.Run("Q", M, M, M, ref WORK, IR + o_work, LDWRKR
3078: , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3079: IWORK = IE + M;
3080: // *
3081: // * Perform bidiagonal QR iteration, computing left
3082: // * singular vectors of L in WORK(IR) and computing
3083: // * right singular vectors of L in WORK(IU)
3084: // * (Workspace: need 2*M*M+BDSPAC)
3085: // *
3086: this._dbdsqr.Run("U", M, M, M, 0, ref S, offset_s
3087: , ref WORK, IE + o_work, ref WORK, IU + o_work, LDWRKU, ref WORK, IR + o_work, LDWRKR, ref DUM, offset_dum
3088: , 1, ref WORK, IWORK + o_work, ref INFO);
3089: // *
3090: // * Multiply right singular vectors of L in WORK(IU) by
3091: // * Q in A, storing result in VT
3092: // * (Workspace: need M*M)
3093: // *
3094: this._dgemm.Run("N", "N", M, N, M, ONE
3095: , WORK, IU + o_work, LDWRKU, A, offset_a, LDA, ZERO, ref VT, offset_vt
3096: , LDVT);
3097: // *
3098: // * Copy left singular vectors of L to A
3099: // * (Workspace: need M*M)
3100: // *
3101: this._dlacpy.Run("F", M, M, WORK, IR + o_work, LDWRKR, ref A, offset_a
3102: , LDA);
3103: // *
3104: }
3105: else
3106: {
3107: // *
3108: // * Insufficient workspace for a fast algorithm
3109: // *
3110: ITAU = 1;
3111: IWORK = ITAU + M;
3112: // *
3113: // * Compute A=L*Q, copying result to VT
3114: // * (Workspace: need 2*M, prefer M+M*NB)
3115: // *
3116: this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
3117: , LWORK - IWORK + 1, ref IERR);
3118: this._dlacpy.Run("U", M, N, A, offset_a, LDA, ref VT, offset_vt
3119: , LDVT);
3120: // *
3121: // * Generate Q in VT
3122: // * (Workspace: need 2*M, prefer M+M*NB)
3123: // *
3124: this._dorglq.Run(M, N, M, ref VT, offset_vt, LDVT, WORK, ITAU + o_work
3125: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3126: IE = ITAU;
3127: ITAUQ = IE + M;
3128: ITAUP = ITAUQ + M;
3129: IWORK = ITAUP + M;
3130: // *
3131: // * Zero out above L in A
3132: // *
3133: this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref A, 1+2 * LDA + o_a
3134: , LDA);
3135: // *
3136: // * Bidiagonalize L in A
3137: // * (Workspace: need 4*M, prefer 3*M+2*M*NB)
3138: // *
3139: this._dgebrd.Run(M, M, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
3140: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3141: // *
3142: // * Multiply right vectors bidiagonalizing L by Q in VT
3143: // * (Workspace: need 3*M+N, prefer 3*M+N*NB)
3144: // *
3145: this._dormbr.Run("P", "L", "T", M, N, M
3146: , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, IWORK + o_work
3147: , LWORK - IWORK + 1, ref IERR);
3148: // *
3149: // * Generate left bidiagonalizing vectors of L in A
3150: // * (Workspace: need 4*M, prefer 3*M+M*NB)
3151: // *
3152: this._dorgbr.Run("Q", M, M, M, ref A, offset_a, LDA
3153: , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3154: IWORK = IE + M;
3155: // *
3156: // * Perform bidiagonal QR iteration, compute left
3157: // * singular vectors of A in A and compute right
3158: // * singular vectors of A in VT
3159: // * (Workspace: need BDSPAC)
3160: // *
3161: this._dbdsqr.Run("U", M, N, M, 0, ref S, offset_s
3162: , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref A, offset_a, LDA, ref DUM, offset_dum
3163: , 1, ref WORK, IWORK + o_work, ref INFO);
3164: // *
3165: }
3166: // *
3167: }
3168: else
3169: {
3170: if (WNTUAS)
3171: {
3172: // *
3173: // * Path 6t(N much larger than M, JOBU='S' or 'A',
3174: // * JOBVT='S')
3175: // * M right singular vectors to be computed in VT and
3176: // * M left singular vectors to be computed in U
3177: // *
3178: if (LWORK >= M * M + Math.Max(4 * M, BDSPAC))
3179: {
3180: // *
3181: // * Sufficient workspace for a fast algorithm
3182: // *
3183: IU = 1;
3184: if (LWORK >= WRKBL + LDA * M)
3185: {
3186: // *
3187: // * WORK(IU) is LDA by N
3188: // *
3189: LDWRKU = LDA;
3190: }
3191: else
3192: {
3193: // *
3194: // * WORK(IU) is LDA by M
3195: // *
3196: LDWRKU = M;
3197: }
3198: ITAU = IU + LDWRKU * M;
3199: IWORK = ITAU + M;
3200: // *
3201: // * Compute A=L*Q
3202: // * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
3203: // *
3204: this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
3205: , LWORK - IWORK + 1, ref IERR);
3206: // *
3207: // * Copy L to WORK(IU), zeroing out above it
3208: // *
3209: this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref WORK, IU + o_work
3210: , LDWRKU);
3211: this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref WORK, IU + LDWRKU + o_work
3212: , LDWRKU);
3213: // *
3214: // * Generate Q in A
3215: // * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
3216: // *
3217: this._dorglq.Run(M, N, M, ref A, offset_a, LDA, WORK, ITAU + o_work
3218: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3219: IE = ITAU;
3220: ITAUQ = IE + M;
3221: ITAUP = ITAUQ + M;
3222: IWORK = ITAUP + M;
3223: // *
3224: // * Bidiagonalize L in WORK(IU), copying result to U
3225: // * (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
3226: // *
3227: this._dgebrd.Run(M, M, ref WORK, IU + o_work, LDWRKU, ref S, offset_s, ref WORK, IE + o_work
3228: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3229: this._dlacpy.Run("L", M, M, WORK, IU + o_work, LDWRKU, ref U, offset_u
3230: , LDU);
3231: // *
3232: // * Generate right bidiagonalizing vectors in WORK(IU)
3233: // * (Workspace: need M*M+4*M-1,
3234: // * prefer M*M+3*M+(M-1)*NB)
3235: // *
3236: this._dorgbr.Run("P", M, M, M, ref WORK, IU + o_work, LDWRKU
3237: , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3238: // *
3239: // * Generate left bidiagonalizing vectors in U
3240: // * (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB)
3241: // *
3242: this._dorgbr.Run("Q", M, M, M, ref U, offset_u, LDU
3243: , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3244: IWORK = IE + M;
3245: // *
3246: // * Perform bidiagonal QR iteration, computing left
3247: // * singular vectors of L in U and computing right
3248: // * singular vectors of L in WORK(IU)
3249: // * (Workspace: need M*M+BDSPAC)
3250: // *
3251: this._dbdsqr.Run("U", M, M, M, 0, ref S, offset_s
3252: , ref WORK, IE + o_work, ref WORK, IU + o_work, LDWRKU, ref U, offset_u, LDU, ref DUM, offset_dum
3253: , 1, ref WORK, IWORK + o_work, ref INFO);
3254: // *
3255: // * Multiply right singular vectors of L in WORK(IU) by
3256: // * Q in A, storing result in VT
3257: // * (Workspace: need M*M)
3258: // *
3259: this._dgemm.Run("N", "N", M, N, M, ONE
3260: , WORK, IU + o_work, LDWRKU, A, offset_a, LDA, ZERO, ref VT, offset_vt
3261: , LDVT);
3262: // *
3263: }
3264: else
3265: {
3266: // *
3267: // * Insufficient workspace for a fast algorithm
3268: // *
3269: ITAU = 1;
3270: IWORK = ITAU + M;
3271: // *
3272: // * Compute A=L*Q, copying result to VT
3273: // * (Workspace: need 2*M, prefer M+M*NB)
3274: // *
3275: this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
3276: , LWORK - IWORK + 1, ref IERR);
3277: this._dlacpy.Run("U", M, N, A, offset_a, LDA, ref VT, offset_vt
3278: , LDVT);
3279: // *
3280: // * Generate Q in VT
3281: // * (Workspace: need 2*M, prefer M+M*NB)
3282: // *
3283: this._dorglq.Run(M, N, M, ref VT, offset_vt, LDVT, WORK, ITAU + o_work
3284: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3285: // *
3286: // * Copy L to U, zeroing out above it
3287: // *
3288: this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref U, offset_u
3289: , LDU);
3290: this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref U, 1+2 * LDU + o_u
3291: , LDU);
3292: IE = ITAU;
3293: ITAUQ = IE + M;
3294: ITAUP = ITAUQ + M;
3295: IWORK = ITAUP + M;
3296: // *
3297: // * Bidiagonalize L in U
3298: // * (Workspace: need 4*M, prefer 3*M+2*M*NB)
3299: // *
3300: this._dgebrd.Run(M, M, ref U, offset_u, LDU, ref S, offset_s, ref WORK, IE + o_work
3301: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3302: // *
3303: // * Multiply right bidiagonalizing vectors in U by Q
3304: // * in VT
3305: // * (Workspace: need 3*M+N, prefer 3*M+N*NB)
3306: // *
3307: this._dormbr.Run("P", "L", "T", M, N, M
3308: , ref U, offset_u, LDU, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, IWORK + o_work
3309: , LWORK - IWORK + 1, ref IERR);
3310: // *
3311: // * Generate left bidiagonalizing vectors in U
3312: // * (Workspace: need 4*M, prefer 3*M+M*NB)
3313: // *
3314: this._dorgbr.Run("Q", M, M, M, ref U, offset_u, LDU
3315: , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3316: IWORK = IE + M;
3317: // *
3318: // * Perform bidiagonal QR iteration, computing left
3319: // * singular vectors of A in U and computing right
3320: // * singular vectors of A in VT
3321: // * (Workspace: need BDSPAC)
3322: // *
3323: this._dbdsqr.Run("U", M, N, M, 0, ref S, offset_s
3324: , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref U, offset_u, LDU, ref DUM, offset_dum
3325: , 1, ref WORK, IWORK + o_work, ref INFO);
3326: // *
3327: }
3328: // *
3329: }
3330: }
3331: }
3332: // *
3333: }
3334: else
3335: {
3336: if (WNTVA)
3337: {
3338: // *
3339: if (WNTUN)
3340: {
3341: // *
3342: // * Path 7t(N much larger than M, JOBU='N', JOBVT='A')
3343: // * N right singular vectors to be computed in VT and
3344: // * no left singular vectors to be computed
3345: // *
3346: if (LWORK >= M * M + Math.Max(N + M, Math.Max(4 * M, BDSPAC)))
3347: {
3348: // *
3349: // * Sufficient workspace for a fast algorithm
3350: // *
3351: IR = 1;
3352: if (LWORK >= WRKBL + LDA * M)
3353: {
3354: // *
3355: // * WORK(IR) is LDA by M
3356: // *
3357: LDWRKR = LDA;
3358: }
3359: else
3360: {
3361: // *
3362: // * WORK(IR) is M by M
3363: // *
3364: LDWRKR = M;
3365: }
3366: ITAU = IR + LDWRKR * M;
3367: IWORK = ITAU + M;
3368: // *
3369: // * Compute A=L*Q, copying result to VT
3370: // * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
3371: // *
3372: this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
3373: , LWORK - IWORK + 1, ref IERR);
3374: this._dlacpy.Run("U", M, N, A, offset_a, LDA, ref VT, offset_vt
3375: , LDVT);
3376: // *
3377: // * Copy L to WORK(IR), zeroing out above it
3378: // *
3379: this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref WORK, IR + o_work
3380: , LDWRKR);
3381: this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref WORK, IR + LDWRKR + o_work
3382: , LDWRKR);
3383: // *
3384: // * Generate Q in VT
3385: // * (Workspace: need M*M+M+N, prefer M*M+M+N*NB)
3386: // *
3387: this._dorglq.Run(N, N, M, ref VT, offset_vt, LDVT, WORK, ITAU + o_work
3388: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3389: IE = ITAU;
3390: ITAUQ = IE + M;
3391: ITAUP = ITAUQ + M;
3392: IWORK = ITAUP + M;
3393: // *
3394: // * Bidiagonalize L in WORK(IR)
3395: // * (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
3396: // *
3397: this._dgebrd.Run(M, M, ref WORK, IR + o_work, LDWRKR, ref S, offset_s, ref WORK, IE + o_work
3398: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3399: // *
3400: // * Generate right bidiagonalizing vectors in WORK(IR)
3401: // * (Workspace: need M*M+4*M-1,
3402: // * prefer M*M+3*M+(M-1)*NB)
3403: // *
3404: this._dorgbr.Run("P", M, M, M, ref WORK, IR + o_work, LDWRKR
3405: , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3406: IWORK = IE + M;
3407: // *
3408: // * Perform bidiagonal QR iteration, computing right
3409: // * singular vectors of L in WORK(IR)
3410: // * (Workspace: need M*M+BDSPAC)
3411: // *
3412: this._dbdsqr.Run("U", M, M, 0, 0, ref S, offset_s
3413: , ref WORK, IE + o_work, ref WORK, IR + o_work, LDWRKR, ref DUM, offset_dum, 1, ref DUM, offset_dum
3414: , 1, ref WORK, IWORK + o_work, ref INFO);
3415: // *
3416: // * Multiply right singular vectors of L in WORK(IR) by
3417: // * Q in VT, storing result in A
3418: // * (Workspace: need M*M)
3419: // *
3420: this._dgemm.Run("N", "N", M, N, M, ONE
3421: , WORK, IR + o_work, LDWRKR, VT, offset_vt, LDVT, ZERO, ref A, offset_a
3422: , LDA);
3423: // *
3424: // * Copy right singular vectors of A from A to VT
3425: // *
3426: this._dlacpy.Run("F", M, N, A, offset_a, LDA, ref VT, offset_vt
3427: , LDVT);
3428: // *
3429: }
3430: else
3431: {
3432: // *
3433: // * Insufficient workspace for a fast algorithm
3434: // *
3435: ITAU = 1;
3436: IWORK = ITAU + M;
3437: // *
3438: // * Compute A=L*Q, copying result to VT
3439: // * (Workspace: need 2*M, prefer M+M*NB)
3440: // *
3441: this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
3442: , LWORK - IWORK + 1, ref IERR);
3443: this._dlacpy.Run("U", M, N, A, offset_a, LDA, ref VT, offset_vt
3444: , LDVT);
3445: // *
3446: // * Generate Q in VT
3447: // * (Workspace: need M+N, prefer M+N*NB)
3448: // *
3449: this._dorglq.Run(N, N, M, ref VT, offset_vt, LDVT, WORK, ITAU + o_work
3450: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3451: IE = ITAU;
3452: ITAUQ = IE + M;
3453: ITAUP = ITAUQ + M;
3454: IWORK = ITAUP + M;
3455: // *
3456: // * Zero out above L in A
3457: // *
3458: this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref A, 1+2 * LDA + o_a
3459: , LDA);
3460: // *
3461: // * Bidiagonalize L in A
3462: // * (Workspace: need 4*M, prefer 3*M+2*M*NB)
3463: // *
3464: this._dgebrd.Run(M, M, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
3465: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3466: // *
3467: // * Multiply right bidiagonalizing vectors in A by Q
3468: // * in VT
3469: // * (Workspace: need 3*M+N, prefer 3*M+N*NB)
3470: // *
3471: this._dormbr.Run("P", "L", "T", M, N, M
3472: , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, IWORK + o_work
3473: , LWORK - IWORK + 1, ref IERR);
3474: IWORK = IE + M;
3475: // *
3476: // * Perform bidiagonal QR iteration, computing right
3477: // * singular vectors of A in VT
3478: // * (Workspace: need BDSPAC)
3479: // *
3480: this._dbdsqr.Run("U", M, N, 0, 0, ref S, offset_s
3481: , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref DUM, offset_dum, 1, ref DUM, offset_dum
3482: , 1, ref WORK, IWORK + o_work, ref INFO);
3483: // *
3484: }
3485: // *
3486: }
3487: else
3488: {
3489: if (WNTUO)
3490: {
3491: // *
3492: // * Path 8t(N much larger than M, JOBU='O', JOBVT='A')
3493: // * N right singular vectors to be computed in VT and
3494: // * M left singular vectors to be overwritten on A
3495: // *
3496: if (LWORK >= 2 * M * M + Math.Max(N + M, Math.Max(4 * M, BDSPAC)))
3497: {
3498: // *
3499: // * Sufficient workspace for a fast algorithm
3500: // *
3501: IU = 1;
3502: if (LWORK >= WRKBL + 2 * LDA * M)
3503: {
3504: // *
3505: // * WORK(IU) is LDA by M and WORK(IR) is LDA by M
3506: // *
3507: LDWRKU = LDA;
3508: IR = IU + LDWRKU * M;
3509: LDWRKR = LDA;
3510: }
3511: else
3512: {
3513: if (LWORK >= WRKBL + (LDA + M) * M)
3514: {
3515: // *
3516: // * WORK(IU) is LDA by M and WORK(IR) is M by M
3517: // *
3518: LDWRKU = LDA;
3519: IR = IU + LDWRKU * M;
3520: LDWRKR = M;
3521: }
3522: else
3523: {
3524: // *
3525: // * WORK(IU) is M by M and WORK(IR) is M by M
3526: // *
3527: LDWRKU = M;
3528: IR = IU + LDWRKU * M;
3529: LDWRKR = M;
3530: }
3531: }
3532: ITAU = IR + LDWRKR * M;
3533: IWORK = ITAU + M;
3534: // *
3535: // * Compute A=L*Q, copying result to VT
3536: // * (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
3537: // *
3538: this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
3539: , LWORK - IWORK + 1, ref IERR);
3540: this._dlacpy.Run("U", M, N, A, offset_a, LDA, ref VT, offset_vt
3541: , LDVT);
3542: // *
3543: // * Generate Q in VT
3544: // * (Workspace: need 2*M*M+M+N, prefer 2*M*M+M+N*NB)
3545: // *
3546: this._dorglq.Run(N, N, M, ref VT, offset_vt, LDVT, WORK, ITAU + o_work
3547: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3548: // *
3549: // * Copy L to WORK(IU), zeroing out above it
3550: // *
3551: this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref WORK, IU + o_work
3552: , LDWRKU);
3553: this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref WORK, IU + LDWRKU + o_work
3554: , LDWRKU);
3555: IE = ITAU;
3556: ITAUQ = IE + M;
3557: ITAUP = ITAUQ + M;
3558: IWORK = ITAUP + M;
3559: // *
3560: // * Bidiagonalize L in WORK(IU), copying result to
3561: // * WORK(IR)
3562: // * (Workspace: need 2*M*M+4*M,
3563: // * prefer 2*M*M+3*M+2*M*NB)
3564: // *
3565: this._dgebrd.Run(M, M, ref WORK, IU + o_work, LDWRKU, ref S, offset_s, ref WORK, IE + o_work
3566: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3567: this._dlacpy.Run("L", M, M, WORK, IU + o_work, LDWRKU, ref WORK, IR + o_work
3568: , LDWRKR);
3569: // *
3570: // * Generate right bidiagonalizing vectors in WORK(IU)
3571: // * (Workspace: need 2*M*M+4*M-1,
3572: // * prefer 2*M*M+3*M+(M-1)*NB)
3573: // *
3574: this._dorgbr.Run("P", M, M, M, ref WORK, IU + o_work, LDWRKU
3575: , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3576: // *
3577: // * Generate left bidiagonalizing vectors in WORK(IR)
3578: // * (Workspace: need 2*M*M+4*M, prefer 2*M*M+3*M+M*NB)
3579: // *
3580: this._dorgbr.Run("Q", M, M, M, ref WORK, IR + o_work, LDWRKR
3581: , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3582: IWORK = IE + M;
3583: // *
3584: // * Perform bidiagonal QR iteration, computing left
3585: // * singular vectors of L in WORK(IR) and computing
3586: // * right singular vectors of L in WORK(IU)
3587: // * (Workspace: need 2*M*M+BDSPAC)
3588: // *
3589: this._dbdsqr.Run("U", M, M, M, 0, ref S, offset_s
3590: , ref WORK, IE + o_work, ref WORK, IU + o_work, LDWRKU, ref WORK, IR + o_work, LDWRKR, ref DUM, offset_dum
3591: , 1, ref WORK, IWORK + o_work, ref INFO);
3592: // *
3593: // * Multiply right singular vectors of L in WORK(IU) by
3594: // * Q in VT, storing result in A
3595: // * (Workspace: need M*M)
3596: // *
3597: this._dgemm.Run("N", "N", M, N, M, ONE
3598: , WORK, IU + o_work, LDWRKU, VT, offset_vt, LDVT, ZERO, ref A, offset_a
3599: , LDA);
3600: // *
3601: // * Copy right singular vectors of A from A to VT
3602: // *
3603: this._dlacpy.Run("F", M, N, A, offset_a, LDA, ref VT, offset_vt
3604: , LDVT);
3605: // *
3606: // * Copy left singular vectors of A from WORK(IR) to A
3607: // *
3608: this._dlacpy.Run("F", M, M, WORK, IR + o_work, LDWRKR, ref A, offset_a
3609: , LDA);
3610: // *
3611: }
3612: else
3613: {
3614: // *
3615: // * Insufficient workspace for a fast algorithm
3616: // *
3617: ITAU = 1;
3618: IWORK = ITAU + M;
3619: // *
3620: // * Compute A=L*Q, copying result to VT
3621: // * (Workspace: need 2*M, prefer M+M*NB)
3622: // *
3623: this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
3624: , LWORK - IWORK + 1, ref IERR);
3625: this._dlacpy.Run("U", M, N, A, offset_a, LDA, ref VT, offset_vt
3626: , LDVT);
3627: // *
3628: // * Generate Q in VT
3629: // * (Workspace: need M+N, prefer M+N*NB)
3630: // *
3631: this._dorglq.Run(N, N, M, ref VT, offset_vt, LDVT, WORK, ITAU + o_work
3632: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3633: IE = ITAU;
3634: ITAUQ = IE + M;
3635: ITAUP = ITAUQ + M;
3636: IWORK = ITAUP + M;
3637: // *
3638: // * Zero out above L in A
3639: // *
3640: this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref A, 1+2 * LDA + o_a
3641: , LDA);
3642: // *
3643: // * Bidiagonalize L in A
3644: // * (Workspace: need 4*M, prefer 3*M+2*M*NB)
3645: // *
3646: this._dgebrd.Run(M, M, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
3647: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3648: // *
3649: // * Multiply right bidiagonalizing vectors in A by Q
3650: // * in VT
3651: // * (Workspace: need 3*M+N, prefer 3*M+N*NB)
3652: // *
3653: this._dormbr.Run("P", "L", "T", M, N, M
3654: , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, IWORK + o_work
3655: , LWORK - IWORK + 1, ref IERR);
3656: // *
3657: // * Generate left bidiagonalizing vectors in A
3658: // * (Workspace: need 4*M, prefer 3*M+M*NB)
3659: // *
3660: this._dorgbr.Run("Q", M, M, M, ref A, offset_a, LDA
3661: , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3662: IWORK = IE + M;
3663: // *
3664: // * Perform bidiagonal QR iteration, computing left
3665: // * singular vectors of A in A and computing right
3666: // * singular vectors of A in VT
3667: // * (Workspace: need BDSPAC)
3668: // *
3669: this._dbdsqr.Run("U", M, N, M, 0, ref S, offset_s
3670: , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref A, offset_a, LDA, ref DUM, offset_dum
3671: , 1, ref WORK, IWORK + o_work, ref INFO);
3672: // *
3673: }
3674: // *
3675: }
3676: else
3677: {
3678: if (WNTUAS)
3679: {
3680: // *
3681: // * Path 9t(N much larger than M, JOBU='S' or 'A',
3682: // * JOBVT='A')
3683: // * N right singular vectors to be computed in VT and
3684: // * M left singular vectors to be computed in U
3685: // *
3686: if (LWORK >= M * M + Math.Max(N + M, Math.Max(4 * M, BDSPAC)))
3687: {
3688: // *
3689: // * Sufficient workspace for a fast algorithm
3690: // *
3691: IU = 1;
3692: if (LWORK >= WRKBL + LDA * M)
3693: {
3694: // *
3695: // * WORK(IU) is LDA by M
3696: // *
3697: LDWRKU = LDA;
3698: }
3699: else
3700: {
3701: // *
3702: // * WORK(IU) is M by M
3703: // *
3704: LDWRKU = M;
3705: }
3706: ITAU = IU + LDWRKU * M;
3707: IWORK = ITAU + M;
3708: // *
3709: // * Compute A=L*Q, copying result to VT
3710: // * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
3711: // *
3712: this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
3713: , LWORK - IWORK + 1, ref IERR);
3714: this._dlacpy.Run("U", M, N, A, offset_a, LDA, ref VT, offset_vt
3715: , LDVT);
3716: // *
3717: // * Generate Q in VT
3718: // * (Workspace: need M*M+M+N, prefer M*M+M+N*NB)
3719: // *
3720: this._dorglq.Run(N, N, M, ref VT, offset_vt, LDVT, WORK, ITAU + o_work
3721: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3722: // *
3723: // * Copy L to WORK(IU), zeroing out above it
3724: // *
3725: this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref WORK, IU + o_work
3726: , LDWRKU);
3727: this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref WORK, IU + LDWRKU + o_work
3728: , LDWRKU);
3729: IE = ITAU;
3730: ITAUQ = IE + M;
3731: ITAUP = ITAUQ + M;
3732: IWORK = ITAUP + M;
3733: // *
3734: // * Bidiagonalize L in WORK(IU), copying result to U
3735: // * (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
3736: // *
3737: this._dgebrd.Run(M, M, ref WORK, IU + o_work, LDWRKU, ref S, offset_s, ref WORK, IE + o_work
3738: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3739: this._dlacpy.Run("L", M, M, WORK, IU + o_work, LDWRKU, ref U, offset_u
3740: , LDU);
3741: // *
3742: // * Generate right bidiagonalizing vectors in WORK(IU)
3743: // * (Workspace: need M*M+4*M, prefer M*M+3*M+(M-1)*NB)
3744: // *
3745: this._dorgbr.Run("P", M, M, M, ref WORK, IU + o_work, LDWRKU
3746: , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3747: // *
3748: // * Generate left bidiagonalizing vectors in U
3749: // * (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB)
3750: // *
3751: this._dorgbr.Run("Q", M, M, M, ref U, offset_u, LDU
3752: , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3753: IWORK = IE + M;
3754: // *
3755: // * Perform bidiagonal QR iteration, computing left
3756: // * singular vectors of L in U and computing right
3757: // * singular vectors of L in WORK(IU)
3758: // * (Workspace: need M*M+BDSPAC)
3759: // *
3760: this._dbdsqr.Run("U", M, M, M, 0, ref S, offset_s
3761: , ref WORK, IE + o_work, ref WORK, IU + o_work, LDWRKU, ref U, offset_u, LDU, ref DUM, offset_dum
3762: , 1, ref WORK, IWORK + o_work, ref INFO);
3763: // *
3764: // * Multiply right singular vectors of L in WORK(IU) by
3765: // * Q in VT, storing result in A
3766: // * (Workspace: need M*M)
3767: // *
3768: this._dgemm.Run("N", "N", M, N, M, ONE
3769: , WORK, IU + o_work, LDWRKU, VT, offset_vt, LDVT, ZERO, ref A, offset_a
3770: , LDA);
3771: // *
3772: // * Copy right singular vectors of A from A to VT
3773: // *
3774: this._dlacpy.Run("F", M, N, A, offset_a, LDA, ref VT, offset_vt
3775: , LDVT);
3776: // *
3777: }
3778: else
3779: {
3780: // *
3781: // * Insufficient workspace for a fast algorithm
3782: // *
3783: ITAU = 1;
3784: IWORK = ITAU + M;
3785: // *
3786: // * Compute A=L*Q, copying result to VT
3787: // * (Workspace: need 2*M, prefer M+M*NB)
3788: // *
3789: this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
3790: , LWORK - IWORK + 1, ref IERR);
3791: this._dlacpy.Run("U", M, N, A, offset_a, LDA, ref VT, offset_vt
3792: , LDVT);
3793: // *
3794: // * Generate Q in VT
3795: // * (Workspace: need M+N, prefer M+N*NB)
3796: // *
3797: this._dorglq.Run(N, N, M, ref VT, offset_vt, LDVT, WORK, ITAU + o_work
3798: , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3799: // *
3800: // * Copy L to U, zeroing out above it
3801: // *
3802: this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref U, offset_u
3803: , LDU);
3804: this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref U, 1+2 * LDU + o_u
3805: , LDU);
3806: IE = ITAU;
3807: ITAUQ = IE + M;
3808: ITAUP = ITAUQ + M;
3809: IWORK = ITAUP + M;
3810: // *
3811: // * Bidiagonalize L in U
3812: // * (Workspace: need 4*M, prefer 3*M+2*M*NB)
3813: // *
3814: this._dgebrd.Run(M, M, ref U, offset_u, LDU, ref S, offset_s, ref WORK, IE + o_work
3815: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3816: // *
3817: // * Multiply right bidiagonalizing vectors in U by Q
3818: // * in VT
3819: // * (Workspace: need 3*M+N, prefer 3*M+N*NB)
3820: // *
3821: this._dormbr.Run("P", "L", "T", M, N, M
3822: , ref U, offset_u, LDU, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, IWORK + o_work
3823: , LWORK - IWORK + 1, ref IERR);
3824: // *
3825: // * Generate left bidiagonalizing vectors in U
3826: // * (Workspace: need 4*M, prefer 3*M+M*NB)
3827: // *
3828: this._dorgbr.Run("Q", M, M, M, ref U, offset_u, LDU
3829: , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3830: IWORK = IE + M;
3831: // *
3832: // * Perform bidiagonal QR iteration, computing left
3833: // * singular vectors of A in U and computing right
3834: // * singular vectors of A in VT
3835: // * (Workspace: need BDSPAC)
3836: // *
3837: this._dbdsqr.Run("U", M, N, M, 0, ref S, offset_s
3838: , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref U, offset_u, LDU, ref DUM, offset_dum
3839: , 1, ref WORK, IWORK + o_work, ref INFO);
3840: // *
3841: }
3842: // *
3843: }
3844: }
3845: }
3846: // *
3847: }
3848: }
3849: }
3850: }
3851: }
3852: // *
3853: }
3854: else
3855: {
3856: // *
3857: // * N .LT. MNTHR
3858: // *
3859: // * Path 10t(N greater than M, but not much larger)
3860: // * Reduce to bidiagonal form without LQ decomposition
3861: // *
3862: IE = 1;
3863: ITAUQ = IE + M;
3864: ITAUP = ITAUQ + M;
3865: IWORK = ITAUP + M;
3866: // *
3867: // * Bidiagonalize A
3868: // * (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
3869: // *
3870: this._dgebrd.Run(M, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
3871: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3872: if (WNTUAS)
3873: {
3874: // *
3875: // * If left singular vectors desired in U, copy result to U
3876: // * and generate left bidiagonalizing vectors in U
3877: // * (Workspace: need 4*M-1, prefer 3*M+(M-1)*NB)
3878: // *
3879: this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref U, offset_u
3880: , LDU);
3881: this._dorgbr.Run("Q", M, M, N, ref U, offset_u, LDU
3882: , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3883: }
3884: if (WNTVAS)
3885: {
3886: // *
3887: // * If right singular vectors desired in VT, copy result to
3888: // * VT and generate right bidiagonalizing vectors in VT
3889: // * (Workspace: need 3*M+NRVT, prefer 3*M+NRVT*NB)
3890: // *
3891: this._dlacpy.Run("U", M, N, A, offset_a, LDA, ref VT, offset_vt
3892: , LDVT);
3893: if (WNTVA) NRVT = N;
3894: if (WNTVS) NRVT = M;
3895: this._dorgbr.Run("P", NRVT, N, M, ref VT, offset_vt, LDVT
3896: , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3897: }
3898: if (WNTUO)
3899: {
3900: // *
3901: // * If left singular vectors desired in A, generate left
3902: // * bidiagonalizing vectors in A
3903: // * (Workspace: need 4*M-1, prefer 3*M+(M-1)*NB)
3904: // *
3905: this._dorgbr.Run("Q", M, M, N, ref A, offset_a, LDA
3906: , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3907: }
3908: if (WNTVO)
3909: {
3910: // *
3911: // * If right singular vectors desired in A, generate right
3912: // * bidiagonalizing vectors in A
3913: // * (Workspace: need 4*M, prefer 3*M+M*NB)
3914: // *
3915: this._dorgbr.Run("P", M, N, M, ref A, offset_a, LDA
3916: , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3917: }
3918: IWORK = IE + M;
3919: if (WNTUAS || WNTUO) NRU = M;
3920: if (WNTUN) NRU = 0;
3921: if (WNTVAS || WNTVO) NCVT = N;
3922: if (WNTVN) NCVT = 0;
3923: if ((!WNTUO) && (!WNTVO))
3924: {
3925: // *
3926: // * Perform bidiagonal QR iteration, if desired, computing
3927: // * left singular vectors in U and computing right singular
3928: // * vectors in VT
3929: // * (Workspace: need BDSPAC)
3930: // *
3931: this._dbdsqr.Run("L", M, NCVT, NRU, 0, ref S, offset_s
3932: , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref U, offset_u, LDU, ref DUM, offset_dum
3933: , 1, ref WORK, IWORK + o_work, ref INFO);
3934: }
3935: else
3936: {
3937: if ((!WNTUO) && WNTVO)
3938: {
3939: // *
3940: // * Perform bidiagonal QR iteration, if desired, computing
3941: // * left singular vectors in U and computing right singular
3942: // * vectors in A
3943: // * (Workspace: need BDSPAC)
3944: // *
3945: this._dbdsqr.Run("L", M, NCVT, NRU, 0, ref S, offset_s
3946: , ref WORK, IE + o_work, ref A, offset_a, LDA, ref U, offset_u, LDU, ref DUM, offset_dum
3947: , 1, ref WORK, IWORK + o_work, ref INFO);
3948: }
3949: else
3950: {
3951: // *
3952: // * Perform bidiagonal QR iteration, if desired, computing
3953: // * left singular vectors in A and computing right singular
3954: // * vectors in VT
3955: // * (Workspace: need BDSPAC)
3956: // *
3957: this._dbdsqr.Run("L", M, NCVT, NRU, 0, ref S, offset_s
3958: , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref A, offset_a, LDA, ref DUM, offset_dum
3959: , 1, ref WORK, IWORK + o_work, ref INFO);
3960: }
3961: }
3962: // *
3963: }
3964: // *
3965: }
3966: // *
3967: // * If DBDSQR failed to converge, copy unconverged superdiagonals
3968: // * to WORK( 2:MINMN )
3969: // *
3970: if (INFO != 0)
3971: {
3972: if (IE > 2)
3973: {
3974: for (I = 1; I <= MINMN - 1; I++)
3975: {
3976: WORK[I + 1 + o_work] = WORK[I + IE - 1 + o_work];
3977: }
3978: }
3979: if (IE < 2)
3980: {
3981: for (I = MINMN - 1; I >= 1; I += - 1)
3982: {
3983: WORK[I + 1 + o_work] = WORK[I + IE - 1 + o_work];
3984: }
3985: }
3986: }
3987: // *
3988: // * Undo scaling if necessary
3989: // *
3990: if (ISCL == 1)
3991: {
3992: if (ANRM > BIGNUM)
3993: {
3994: this._dlascl.Run("G", 0, 0, BIGNUM, ANRM, MINMN
3995: , 1, ref S, offset_s, MINMN, ref IERR);
3996: }
3997: if (INFO != 0 && ANRM > BIGNUM)
3998: {
3999: this._dlascl.Run("G", 0, 0, BIGNUM, ANRM, MINMN - 1
4000: , 1, ref WORK, 2 + o_work, MINMN, ref IERR);
4001: }
4002: if (ANRM < SMLNUM)
4003: {
4004: this._dlascl.Run("G", 0, 0, SMLNUM, ANRM, MINMN
4005: , 1, ref S, offset_s, MINMN, ref IERR);
4006: }
4007: if (INFO != 0 && ANRM < SMLNUM)
4008: {
4009: this._dlascl.Run("G", 0, 0, SMLNUM, ANRM, MINMN - 1
4010: , 1, ref WORK, 2 + o_work, MINMN, ref IERR);
4011: }
4012: }
4013: // *
4014: // * Return optimal workspace in WORK(1)
4015: // *
4016: WORK[1 + o_work] = MAXWRK;
4017: // *
4018: return;
4019: // *
4020: // * End of DGESVD
4021: // *
4022:
4023: #endregion
4024:
4025: }
4026: }
4027: }