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: /// DGELSD computes the minimum-norm solution to a real linear least
27: /// squares problem:
28: /// minimize 2-norm(| b - A*x |)
29: /// using the singular value decomposition (SVD) of A. A is an M-by-N
30: /// matrix which may be rank-deficient.
31: ///
32: /// Several right hand side vectors b and solution vectors x can be
33: /// handled in a single call; they are stored as the columns of the
34: /// M-by-NRHS right hand side matrix B and the N-by-NRHS solution
35: /// matrix X.
36: ///
37: /// The problem is solved in three steps:
38: /// (1) Reduce the coefficient matrix A to bidiagonal form with
39: /// Householder transformations, reducing the original problem
40: /// into a "bidiagonal least squares problem" (BLS)
41: /// (2) Solve the BLS using a divide and conquer approach.
42: /// (3) Apply back all the Householder tranformations to solve
43: /// the original least squares problem.
44: ///
45: /// The effective rank of A is determined by treating as zero those
46: /// singular values which are less than RCOND times the largest singular
47: /// value.
48: ///
49: /// The divide and conquer algorithm makes very mild assumptions about
50: /// floating point arithmetic. It will work on machines with a guard
51: /// digit in add/subtract, or on those binary machines without guard
52: /// digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
53: /// Cray-2. It could conceivably fail on hexadecimal or decimal machines
54: /// without guard digits, but we know of none.
55: ///
56: ///</summary>
57: public class DGELSD
58: {
59:
60:
61: #region Dependencies
62:
63: DGEBRD _dgebrd; DGELQF _dgelqf; DGEQRF _dgeqrf; DLABAD _dlabad; DLACPY _dlacpy; DLALSD _dlalsd; DLASCL _dlascl;
64: DLASET _dlaset;DORMBR _dormbr; DORMLQ _dormlq; DORMQR _dormqr; XERBLA _xerbla; ILAENV _ilaenv; DLAMCH _dlamch;
65: DLANGE _dlange;
66:
67: #endregion
68:
69:
70: #region Fields
71:
72: const double ZERO = 0.0E0; const double ONE = 1.0E0; const double TWO = 2.0E0; bool LQUERY = false; int IASCL = 0;
73: int IBSCL = 0;int IE = 0; int IL = 0; int ITAU = 0; int ITAUP = 0; int ITAUQ = 0; int LDWORK = 0; int MAXMN = 0;
74: int MAXWRK = 0;int MINMN = 0; int MINWRK = 0; int MM = 0; int MNTHR = 0; int NLVL = 0; int NWORK = 0; int SMLSIZ = 0;
75: int WLALSD = 0;double ANRM = 0; double BIGNUM = 0; double BNRM = 0; double EPS = 0; double SFMIN = 0; double SMLNUM = 0;
76:
77: #endregion
78:
79: public DGELSD(DGEBRD dgebrd, DGELQF dgelqf, DGEQRF dgeqrf, DLABAD dlabad, DLACPY dlacpy, DLALSD dlalsd, DLASCL dlascl, DLASET dlaset, DORMBR dormbr, DORMLQ dormlq
80: , DORMQR dormqr, XERBLA xerbla, ILAENV ilaenv, DLAMCH dlamch, DLANGE dlange)
81: {
82:
83:
84: #region Set Dependencies
85:
86: this._dgebrd = dgebrd; this._dgelqf = dgelqf; this._dgeqrf = dgeqrf; this._dlabad = dlabad; this._dlacpy = dlacpy;
87: this._dlalsd = dlalsd;this._dlascl = dlascl; this._dlaset = dlaset; this._dormbr = dormbr; this._dormlq = dormlq;
88: this._dormqr = dormqr;this._xerbla = xerbla; this._ilaenv = ilaenv; this._dlamch = dlamch; this._dlange = dlange;
89:
90: #endregion
91:
92: }
93:
94: public DGELSD()
95: {
96:
97:
98: #region Dependencies (Initialization)
99:
100: LSAME lsame = new LSAME();
101: XERBLA xerbla = new XERBLA();
102: DLAMC3 dlamc3 = new DLAMC3();
103: DLAPY2 dlapy2 = new DLAPY2();
104: DNRM2 dnrm2 = new DNRM2();
105: DSCAL dscal = new DSCAL();
106: IEEECK ieeeck = new IEEECK();
107: IPARMQ iparmq = new IPARMQ();
108: DCOPY dcopy = new DCOPY();
109: DLABAD dlabad = new DLABAD();
110: IDAMAX idamax = new IDAMAX();
111: DLASSQ dlassq = new DLASSQ();
112: DROT drot = new DROT();
113: DLASDT dlasdt = new DLASDT();
114: DLAMRG dlamrg = new DLAMRG();
115: DLASD5 dlasd5 = new DLASD5();
116: DDOT ddot = new DDOT();
117: DLAS2 dlas2 = new DLAS2();
118: DLASQ5 dlasq5 = new DLASQ5();
119: DLAZQ4 dlazq4 = new DLAZQ4();
120: DSWAP dswap = new DSWAP();
121: DGEMV dgemv = new DGEMV(lsame, xerbla);
122: DGER dger = new DGER(xerbla);
123: DLARF dlarf = new DLARF(dgemv, dger, lsame);
124: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
125: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
126: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
127: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
128: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
129: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
130: DGEBD2 dgebd2 = new DGEBD2(dlarf, dlarfg, xerbla);
131: DGEMM dgemm = new DGEMM(lsame, xerbla);
132: DLABRD dlabrd = new DLABRD(dgemv, dlarfg, dscal);
133: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
134: DGEBRD dgebrd = new DGEBRD(dgebd2, dgemm, dlabrd, xerbla, ilaenv);
135: DGELQ2 dgelq2 = new DGELQ2(dlarf, dlarfg, xerbla);
136: DTRMM dtrmm = new DTRMM(lsame, xerbla);
137: DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
138: DTRMV dtrmv = new DTRMV(lsame, xerbla);
139: DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);
140: DGELQF dgelqf = new DGELQF(dgelq2, dlarfb, dlarft, xerbla, ilaenv);
141: DGEQR2 dgeqr2 = new DGEQR2(dlarf, dlarfg, xerbla);
142: DGEQRF dgeqrf = new DGEQRF(dgeqr2, dlarfb, dlarft, xerbla, ilaenv);
143: DLACPY dlacpy = new DLACPY(lsame);
144: DLANST dlanst = new DLANST(lsame, dlassq);
145: DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
146: DLALS0 dlals0 = new DLALS0(dcopy, dgemv, dlacpy, dlascl, drot, dscal, xerbla, dlamc3, dnrm2);
147: DLALSA dlalsa = new DLALSA(dcopy, dgemm, dlals0, dlasdt, xerbla);
148: DLARTG dlartg = new DLARTG(dlamch);
149: DLASD7 dlasd7 = new DLASD7(dcopy, dlamrg, drot, xerbla, dlamch, dlapy2);
150: DLAED6 dlaed6 = new DLAED6(dlamch);
151: DLASD4 dlasd4 = new DLASD4(dlaed6, dlasd5, dlamch);
152: DLASET dlaset = new DLASET(lsame);
153: DLASD8 dlasd8 = new DLASD8(dcopy, dlascl, dlasd4, dlaset, xerbla, ddot, dlamc3, dnrm2);
154: DLASD6 dlasd6 = new DLASD6(dcopy, dlamrg, dlascl, dlasd7, dlasd8, xerbla);
155: DLASQ6 dlasq6 = new DLASQ6(dlamch);
156: DLAZQ3 dlazq3 = new DLAZQ3(dlasq5, dlasq6, dlazq4, dlamch);
157: DLASRT dlasrt = new DLASRT(lsame, xerbla);
158: DLASQ2 dlasq2 = new DLASQ2(dlazq3, dlasrt, xerbla, dlamch, ilaenv);
159: DLASQ1 dlasq1 = new DLASQ1(dcopy, dlas2, dlascl, dlasq2, dlasrt, xerbla, dlamch);
160: DLASR dlasr = new DLASR(lsame, xerbla);
161: DLASV2 dlasv2 = new DLASV2(dlamch);
162: DBDSQR dbdsqr = new DBDSQR(lsame, dlamch, dlartg, dlas2, dlasq1, dlasr, dlasv2, drot, dscal, dswap
163: , xerbla);
164: DLASDQ dlasdq = new DLASDQ(dbdsqr, dlartg, dlasr, dswap, xerbla, lsame);
165: DLASDA dlasda = new DLASDA(dcopy, dlasd6, dlasdq, dlasdt, dlaset, xerbla);
166: DLALSD dlalsd = new DLALSD(idamax, dlamch, dlanst, dcopy, dgemm, dlacpy, dlalsa, dlartg, dlascl, dlasda
167: , dlasdq, dlaset, dlasrt, drot, xerbla);
168: DORML2 dorml2 = new DORML2(lsame, dlarf, xerbla);
169: DORMLQ dormlq = new DORMLQ(lsame, ilaenv, dlarfb, dlarft, dorml2, xerbla);
170: DORM2R dorm2r = new DORM2R(lsame, dlarf, xerbla);
171: DORMQR dormqr = new DORMQR(lsame, ilaenv, dlarfb, dlarft, dorm2r, xerbla);
172: DORMBR dormbr = new DORMBR(lsame, ilaenv, dormlq, dormqr, xerbla);
173: DLANGE dlange = new DLANGE(dlassq, lsame);
174:
175: #endregion
176:
177:
178: #region Set Dependencies
179:
180: this._dgebrd = dgebrd; this._dgelqf = dgelqf; this._dgeqrf = dgeqrf; this._dlabad = dlabad; this._dlacpy = dlacpy;
181: this._dlalsd = dlalsd;this._dlascl = dlascl; this._dlaset = dlaset; this._dormbr = dormbr; this._dormlq = dormlq;
182: this._dormqr = dormqr;this._xerbla = xerbla; this._ilaenv = ilaenv; this._dlamch = dlamch; this._dlange = dlange;
183:
184: #endregion
185:
186: }
187: /// <summary>
188: /// Purpose
189: /// =======
190: ///
191: /// DGELSD computes the minimum-norm solution to a real linear least
192: /// squares problem:
193: /// minimize 2-norm(| b - A*x |)
194: /// using the singular value decomposition (SVD) of A. A is an M-by-N
195: /// matrix which may be rank-deficient.
196: ///
197: /// Several right hand side vectors b and solution vectors x can be
198: /// handled in a single call; they are stored as the columns of the
199: /// M-by-NRHS right hand side matrix B and the N-by-NRHS solution
200: /// matrix X.
201: ///
202: /// The problem is solved in three steps:
203: /// (1) Reduce the coefficient matrix A to bidiagonal form with
204: /// Householder transformations, reducing the original problem
205: /// into a "bidiagonal least squares problem" (BLS)
206: /// (2) Solve the BLS using a divide and conquer approach.
207: /// (3) Apply back all the Householder tranformations to solve
208: /// the original least squares problem.
209: ///
210: /// The effective rank of A is determined by treating as zero those
211: /// singular values which are less than RCOND times the largest singular
212: /// value.
213: ///
214: /// The divide and conquer algorithm makes very mild assumptions about
215: /// floating point arithmetic. It will work on machines with a guard
216: /// digit in add/subtract, or on those binary machines without guard
217: /// digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
218: /// Cray-2. It could conceivably fail on hexadecimal or decimal machines
219: /// without guard digits, but we know of none.
220: ///
221: ///</summary>
222: /// <param name="M">
223: /// (input) INTEGER
224: /// The number of rows of A. M .GE. 0.
225: ///</param>
226: /// <param name="N">
227: /// (input) INTEGER
228: /// The number of columns of A. N .GE. 0.
229: ///</param>
230: /// <param name="NRHS">
231: /// (input) INTEGER
232: /// The number of right hand sides, i.e., the number of columns
233: /// of the matrices B and X. NRHS .GE. 0.
234: ///</param>
235: /// <param name="A">
236: /// (input) DOUBLE PRECISION array, dimension (LDA,N)
237: /// On entry, the M-by-N matrix A.
238: /// On exit, A has been destroyed.
239: ///</param>
240: /// <param name="LDA">
241: /// (input) INTEGER
242: /// The leading dimension of the array A. LDA .GE. max(1,M).
243: ///</param>
244: /// <param name="B">
245: /// (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
246: /// On entry, the M-by-NRHS right hand side matrix B.
247: /// On exit, B is overwritten by the N-by-NRHS solution
248: /// matrix X. If m .GE. n and RANK = n, the residual
249: /// sum-of-squares for the solution in the i-th column is given
250: /// by the sum of squares of elements n+1:m in that column.
251: ///</param>
252: /// <param name="LDB">
253: /// (input) INTEGER
254: /// The leading dimension of the array B. LDB .GE. max(1,max(M,N)).
255: ///</param>
256: /// <param name="S">
257: /// (output) DOUBLE PRECISION array, dimension (min(M,N))
258: /// The singular values of A in decreasing order.
259: /// The condition number of A in the 2-norm = S(1)/S(min(m,n)).
260: ///</param>
261: /// <param name="RCOND">
262: /// (input) DOUBLE PRECISION
263: /// RCOND is used to determine the effective rank of A.
264: /// Singular values S(i) .LE. RCOND*S(1) are treated as zero.
265: /// If RCOND .LT. 0, machine precision is used instead.
266: ///</param>
267: /// <param name="RANK">
268: /// (output) INTEGER
269: /// The effective rank of A, i.e., the number of singular values
270: /// which are greater than RCOND*S(1).
271: ///</param>
272: /// <param name="WORK">
273: /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
274: /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
275: ///</param>
276: /// <param name="LWORK">
277: /// (input) INTEGER
278: /// The dimension of the array WORK. LWORK must be at least 1.
279: /// The exact minimum amount of workspace needed depends on M,
280: /// N and NRHS. As long as LWORK is at least
281: /// 12*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2,
282: /// if M is greater than or equal to N or
283: /// 12*M + 2*M*SMLSIZ + 8*M*NLVL + M*NRHS + (SMLSIZ+1)**2,
284: /// if M is less than N, the code will execute correctly.
285: /// SMLSIZ is returned by ILAENV and is equal to the maximum
286: /// size of the subproblems at the bottom of the computation
287: /// tree (usually about 25), and
288: /// NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
289: /// For good performance, LWORK should generally be larger.
290: ///
291: /// If LWORK = -1, then a workspace query is assumed; the routine
292: /// only calculates the optimal size of the WORK array, returns
293: /// this value as the first entry of the WORK array, and no error
294: /// message related to LWORK is issued by XERBLA.
295: ///</param>
296: /// <param name="IWORK">
297: /// (workspace) INTEGER array, dimension (MAX(1,LIWORK))
298: /// LIWORK .GE. 3 * MINMN * NLVL + 11 * MINMN,
299: /// where MINMN = MIN( M,N ).
300: ///</param>
301: /// <param name="INFO">
302: /// (output) INTEGER
303: /// = 0: successful exit
304: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
305: /// .GT. 0: the algorithm for computing the SVD failed to converge;
306: /// if INFO = i, i off-diagonal elements of an intermediate
307: /// bidiagonal form did not converge to zero.
308: ///</param>
309: public void Run(int M, int N, int NRHS, ref double[] A, int offset_a, int LDA, ref double[] B, int offset_b
310: , int LDB, ref double[] S, int offset_s, double RCOND, ref int RANK, ref double[] WORK, int offset_work, int LWORK
311: , ref int[] IWORK, int offset_iwork, ref int INFO)
312: {
313:
314: #region Array Index Correction
315:
316: int o_a = -1 - LDA + offset_a; int o_b = -1 - LDB + offset_b; int o_s = -1 + offset_s;
317: int o_work = -1 + offset_work; int o_iwork = -1 + offset_iwork;
318:
319: #endregion
320:
321:
322: #region Prolog
323:
324: // *
325: // * -- LAPACK driver routine (version 3.1) --
326: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
327: // * November 2006
328: // *
329: // * .. Scalar Arguments ..
330: // * ..
331: // * .. Array Arguments ..
332: // * ..
333: // *
334: // * Purpose
335: // * =======
336: // *
337: // * DGELSD computes the minimum-norm solution to a real linear least
338: // * squares problem:
339: // * minimize 2-norm(| b - A*x |)
340: // * using the singular value decomposition (SVD) of A. A is an M-by-N
341: // * matrix which may be rank-deficient.
342: // *
343: // * Several right hand side vectors b and solution vectors x can be
344: // * handled in a single call; they are stored as the columns of the
345: // * M-by-NRHS right hand side matrix B and the N-by-NRHS solution
346: // * matrix X.
347: // *
348: // * The problem is solved in three steps:
349: // * (1) Reduce the coefficient matrix A to bidiagonal form with
350: // * Householder transformations, reducing the original problem
351: // * into a "bidiagonal least squares problem" (BLS)
352: // * (2) Solve the BLS using a divide and conquer approach.
353: // * (3) Apply back all the Householder tranformations to solve
354: // * the original least squares problem.
355: // *
356: // * The effective rank of A is determined by treating as zero those
357: // * singular values which are less than RCOND times the largest singular
358: // * value.
359: // *
360: // * The divide and conquer algorithm makes very mild assumptions about
361: // * floating point arithmetic. It will work on machines with a guard
362: // * digit in add/subtract, or on those binary machines without guard
363: // * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
364: // * Cray-2. It could conceivably fail on hexadecimal or decimal machines
365: // * without guard digits, but we know of none.
366: // *
367: // * Arguments
368: // * =========
369: // *
370: // * M (input) INTEGER
371: // * The number of rows of A. M >= 0.
372: // *
373: // * N (input) INTEGER
374: // * The number of columns of A. N >= 0.
375: // *
376: // * NRHS (input) INTEGER
377: // * The number of right hand sides, i.e., the number of columns
378: // * of the matrices B and X. NRHS >= 0.
379: // *
380: // * A (input) DOUBLE PRECISION array, dimension (LDA,N)
381: // * On entry, the M-by-N matrix A.
382: // * On exit, A has been destroyed.
383: // *
384: // * LDA (input) INTEGER
385: // * The leading dimension of the array A. LDA >= max(1,M).
386: // *
387: // * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
388: // * On entry, the M-by-NRHS right hand side matrix B.
389: // * On exit, B is overwritten by the N-by-NRHS solution
390: // * matrix X. If m >= n and RANK = n, the residual
391: // * sum-of-squares for the solution in the i-th column is given
392: // * by the sum of squares of elements n+1:m in that column.
393: // *
394: // * LDB (input) INTEGER
395: // * The leading dimension of the array B. LDB >= max(1,max(M,N)).
396: // *
397: // * S (output) DOUBLE PRECISION array, dimension (min(M,N))
398: // * The singular values of A in decreasing order.
399: // * The condition number of A in the 2-norm = S(1)/S(min(m,n)).
400: // *
401: // * RCOND (input) DOUBLE PRECISION
402: // * RCOND is used to determine the effective rank of A.
403: // * Singular values S(i) <= RCOND*S(1) are treated as zero.
404: // * If RCOND < 0, machine precision is used instead.
405: // *
406: // * RANK (output) INTEGER
407: // * The effective rank of A, i.e., the number of singular values
408: // * which are greater than RCOND*S(1).
409: // *
410: // * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
411: // * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
412: // *
413: // * LWORK (input) INTEGER
414: // * The dimension of the array WORK. LWORK must be at least 1.
415: // * The exact minimum amount of workspace needed depends on M,
416: // * N and NRHS. As long as LWORK is at least
417: // * 12*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2,
418: // * if M is greater than or equal to N or
419: // * 12*M + 2*M*SMLSIZ + 8*M*NLVL + M*NRHS + (SMLSIZ+1)**2,
420: // * if M is less than N, the code will execute correctly.
421: // * SMLSIZ is returned by ILAENV and is equal to the maximum
422: // * size of the subproblems at the bottom of the computation
423: // * tree (usually about 25), and
424: // * NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
425: // * For good performance, LWORK should generally be larger.
426: // *
427: // * If LWORK = -1, then a workspace query is assumed; the routine
428: // * only calculates the optimal size of the WORK array, returns
429: // * this value as the first entry of the WORK array, and no error
430: // * message related to LWORK is issued by XERBLA.
431: // *
432: // * IWORK (workspace) INTEGER array, dimension (MAX(1,LIWORK))
433: // * LIWORK >= 3 * MINMN * NLVL + 11 * MINMN,
434: // * where MINMN = MIN( M,N ).
435: // *
436: // * INFO (output) INTEGER
437: // * = 0: successful exit
438: // * < 0: if INFO = -i, the i-th argument had an illegal value.
439: // * > 0: the algorithm for computing the SVD failed to converge;
440: // * if INFO = i, i off-diagonal elements of an intermediate
441: // * bidiagonal form did not converge to zero.
442: // *
443: // * Further Details
444: // * ===============
445: // *
446: // * Based on contributions by
447: // * Ming Gu and Ren-Cang Li, Computer Science Division, University of
448: // * California at Berkeley, USA
449: // * Osni Marques, LBNL/NERSC, USA
450: // *
451: // * =====================================================================
452: // *
453: // * .. Parameters ..
454: // * ..
455: // * .. Local Scalars ..
456: // * ..
457: // * .. External Subroutines ..
458: // * ..
459: // * .. External Functions ..
460: // * ..
461: // * .. Intrinsic Functions ..
462: // INTRINSIC DBLE, INT, LOG, MAX, MIN;
463: // * ..
464: // * .. Executable Statements ..
465: // *
466: // * Test the input arguments.
467: // *
468:
469: #endregion
470:
471:
472: #region Body
473:
474: INFO = 0;
475: MINMN = Math.Min(M, N);
476: MAXMN = Math.Max(M, N);
477: MNTHR = this._ilaenv.Run(6, "DGELSD", " ", M, N, NRHS, - 1);
478: LQUERY = (LWORK == - 1);
479: if (M < 0)
480: {
481: INFO = - 1;
482: }
483: else
484: {
485: if (N < 0)
486: {
487: INFO = - 2;
488: }
489: else
490: {
491: if (NRHS < 0)
492: {
493: INFO = - 3;
494: }
495: else
496: {
497: if (LDA < Math.Max(1, M))
498: {
499: INFO = - 5;
500: }
501: else
502: {
503: if (LDB < Math.Max(1, MAXMN))
504: {
505: INFO = - 7;
506: }
507: }
508: }
509: }
510: }
511: // *
512: SMLSIZ = this._ilaenv.Run(9, "DGELSD", " ", 0, 0, 0, 0);
513: // *
514: // * Compute workspace.
515: // * (Note: Comments in the code beginning "Workspace:" describe the
516: // * minimal amount of workspace needed at that point in the code,
517: // * as well as the preferred amount for good performance.
518: // * NB refers to the optimal block size for the immediately
519: // * following subroutine, as returned by ILAENV.)
520: // *
521: MINWRK = 1;
522: MINMN = Math.Max(1, MINMN);
523: NLVL = Math.Max(Convert.ToInt32(Math.Truncate(Math.Log(Convert.ToDouble(MINMN) / Convert.ToDouble(SMLSIZ + 1)) / Math.Log(TWO))) + 1, 0);
524: // *
525: if (INFO == 0)
526: {
527: MAXWRK = 0;
528: MM = M;
529: if (M >= N && M >= MNTHR)
530: {
531: // *
532: // * Path 1a - overdetermined, with many more rows than columns.
533: // *
534: MM = N;
535: MAXWRK = Math.Max(MAXWRK, N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N, - 1, - 1));
536: MAXWRK = Math.Max(MAXWRK, N + NRHS * this._ilaenv.Run(1, "DORMQR", "LT", M, NRHS, N, - 1));
537: }
538: if (M >= N)
539: {
540: // *
541: // * Path 1 - overdetermined or exactly determined.
542: // *
543: MAXWRK = Math.Max(MAXWRK, 3 * N + (MM + N) * this._ilaenv.Run(1, "DGEBRD", " ", MM, N, - 1, - 1));
544: MAXWRK = Math.Max(MAXWRK, 3 * N + NRHS * this._ilaenv.Run(1, "DORMBR", "QLT", MM, NRHS, N, - 1));
545: MAXWRK = Math.Max(MAXWRK, 3 * N + (N - 1) * this._ilaenv.Run(1, "DORMBR", "PLN", N, NRHS, N, - 1));
546: WLALSD = 9 * N + 2 * N * SMLSIZ + 8 * N * NLVL + N * NRHS + (int)Math.Pow(SMLSIZ + 1, 2);
547: MAXWRK = Math.Max(MAXWRK, 3 * N + WLALSD);
548: MINWRK = Math.Max(3 * N + MM, Math.Max(3 * N + NRHS, 3 * N + WLALSD));
549: }
550: if (N > M)
551: {
552: WLALSD = 9 * M + 2 * M * SMLSIZ + 8 * M * NLVL + M * NRHS + (int)Math.Pow(SMLSIZ + 1, 2);
553: if (N >= MNTHR)
554: {
555: // *
556: // * Path 2a - underdetermined, with many more columns
557: // * than rows.
558: // *
559: MAXWRK = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N, - 1, - 1);
560: MAXWRK = Math.Max(MAXWRK, M * M + 4 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M, - 1, - 1));
561: MAXWRK = Math.Max(MAXWRK, M * M + 4 * M + NRHS * this._ilaenv.Run(1, "DORMBR", "QLT", M, NRHS, M, - 1));
562: MAXWRK = Math.Max(MAXWRK, M * M + 4 * M + (M - 1) * this._ilaenv.Run(1, "DORMBR", "PLN", M, NRHS, M, - 1));
563: if (NRHS > 1)
564: {
565: MAXWRK = Math.Max(MAXWRK, M * M + M + M * NRHS);
566: }
567: else
568: {
569: MAXWRK = Math.Max(MAXWRK, M * M + 2 * M);
570: }
571: MAXWRK = Math.Max(MAXWRK, M + NRHS * this._ilaenv.Run(1, "DORMLQ", "LT", N, NRHS, M, - 1));
572: MAXWRK = Math.Max(MAXWRK, M * M + 4 * M + WLALSD);
573: }
574: else
575: {
576: // *
577: // * Path 2 - remaining underdetermined cases.
578: // *
579: MAXWRK = 3 * M + (N + M) * this._ilaenv.Run(1, "DGEBRD", " ", M, N, - 1, - 1);
580: MAXWRK = Math.Max(MAXWRK, 3 * M + NRHS * this._ilaenv.Run(1, "DORMBR", "QLT", M, NRHS, N, - 1));
581: MAXWRK = Math.Max(MAXWRK, 3 * M + M * this._ilaenv.Run(1, "DORMBR", "PLN", N, NRHS, M, - 1));
582: MAXWRK = Math.Max(MAXWRK, 3 * M + WLALSD);
583: }
584: MINWRK = Math.Max(3 * M + NRHS, Math.Max(3 * M + M, 3 * M + WLALSD));
585: }
586: MINWRK = Math.Min(MINWRK, MAXWRK);
587: WORK[1 + o_work] = MAXWRK;
588: if (LWORK < MINWRK && !LQUERY)
589: {
590: INFO = - 12;
591: }
592: }
593: // *
594: if (INFO != 0)
595: {
596: this._xerbla.Run("DGELSD", - INFO);
597: return;
598: }
599: else
600: {
601: if (LQUERY)
602: {
603: goto LABEL10;
604: }
605: }
606: // *
607: // * Quick return if possible.
608: // *
609: if (M == 0 || N == 0)
610: {
611: RANK = 0;
612: return;
613: }
614: // *
615: // * Get machine parameters.
616: // *
617: EPS = this._dlamch.Run("P");
618: SFMIN = this._dlamch.Run("S");
619: SMLNUM = SFMIN / EPS;
620: BIGNUM = ONE / SMLNUM;
621: this._dlabad.Run(ref SMLNUM, ref BIGNUM);
622: // *
623: // * Scale A if max entry outside range [SMLNUM,BIGNUM].
624: // *
625: ANRM = this._dlange.Run("M", M, N, A, offset_a, LDA, ref WORK, offset_work);
626: IASCL = 0;
627: if (ANRM > ZERO && ANRM < SMLNUM)
628: {
629: // *
630: // * Scale matrix norm up to SMLNUM.
631: // *
632: this._dlascl.Run("G", 0, 0, ANRM, SMLNUM, M
633: , N, ref A, offset_a, LDA, ref INFO);
634: IASCL = 1;
635: }
636: else
637: {
638: if (ANRM > BIGNUM)
639: {
640: // *
641: // * Scale matrix norm down to BIGNUM.
642: // *
643: this._dlascl.Run("G", 0, 0, ANRM, BIGNUM, M
644: , N, ref A, offset_a, LDA, ref INFO);
645: IASCL = 2;
646: }
647: else
648: {
649: if (ANRM == ZERO)
650: {
651: // *
652: // * Matrix all zero. Return zero solution.
653: // *
654: this._dlaset.Run("F", Math.Max(M, N), NRHS, ZERO, ZERO, ref B, offset_b
655: , LDB);
656: this._dlaset.Run("F", MINMN, 1, ZERO, ZERO, ref S, offset_s
657: , 1);
658: RANK = 0;
659: goto LABEL10;
660: }
661: }
662: }
663: // *
664: // * Scale B if max entry outside range [SMLNUM,BIGNUM].
665: // *
666: BNRM = this._dlange.Run("M", M, NRHS, B, offset_b, LDB, ref WORK, offset_work);
667: IBSCL = 0;
668: if (BNRM > ZERO && BNRM < SMLNUM)
669: {
670: // *
671: // * Scale matrix norm up to SMLNUM.
672: // *
673: this._dlascl.Run("G", 0, 0, BNRM, SMLNUM, M
674: , NRHS, ref B, offset_b, LDB, ref INFO);
675: IBSCL = 1;
676: }
677: else
678: {
679: if (BNRM > BIGNUM)
680: {
681: // *
682: // * Scale matrix norm down to BIGNUM.
683: // *
684: this._dlascl.Run("G", 0, 0, BNRM, BIGNUM, M
685: , NRHS, ref B, offset_b, LDB, ref INFO);
686: IBSCL = 2;
687: }
688: }
689: // *
690: // * If M < N make sure certain entries of B are zero.
691: // *
692: if (M < N)
693: {
694: this._dlaset.Run("F", N - M, NRHS, ZERO, ZERO, ref B, M + 1+1 * LDB + o_b
695: , LDB);
696: }
697: // *
698: // * Overdetermined case.
699: // *
700: if (M >= N)
701: {
702: // *
703: // * Path 1 - overdetermined or exactly determined.
704: // *
705: MM = M;
706: if (M >= MNTHR)
707: {
708: // *
709: // * Path 1a - overdetermined, with many more rows than columns.
710: // *
711: MM = N;
712: ITAU = 1;
713: NWORK = ITAU + N;
714: // *
715: // * Compute A=Q*R.
716: // * (Workspace: need 2*N, prefer N+N*NB)
717: // *
718: this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, NWORK + o_work
719: , LWORK - NWORK + 1, ref INFO);
720: // *
721: // * Multiply B by transpose(Q).
722: // * (Workspace: need N+NRHS, prefer N+NRHS*NB)
723: // *
724: this._dormqr.Run("L", "T", M, NRHS, N, ref A, offset_a
725: , LDA, WORK, ITAU + o_work, ref B, offset_b, LDB, ref WORK, NWORK + o_work, LWORK - NWORK + 1
726: , ref INFO);
727: // *
728: // * Zero out below R.
729: // *
730: if (N > 1)
731: {
732: this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref A, 2+1 * LDA + o_a
733: , LDA);
734: }
735: }
736: // *
737: IE = 1;
738: ITAUQ = IE + N;
739: ITAUP = ITAUQ + N;
740: NWORK = ITAUP + N;
741: // *
742: // * Bidiagonalize R in A.
743: // * (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB)
744: // *
745: this._dgebrd.Run(MM, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
746: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref INFO);
747: // *
748: // * Multiply B by transpose of left bidiagonalizing vectors of R.
749: // * (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB)
750: // *
751: this._dormbr.Run("Q", "L", "T", MM, NRHS, N
752: , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref B, offset_b, LDB, ref WORK, NWORK + o_work
753: , LWORK - NWORK + 1, ref INFO);
754: // *
755: // * Solve the bidiagonal least squares problem.
756: // *
757: this._dlalsd.Run("U", SMLSIZ, N, NRHS, ref S, offset_s, ref WORK, IE + o_work
758: , ref B, offset_b, LDB, RCOND, ref RANK, ref WORK, NWORK + o_work, ref IWORK, offset_iwork
759: , ref INFO);
760: if (INFO != 0)
761: {
762: goto LABEL10;
763: }
764: // *
765: // * Multiply B by right bidiagonalizing vectors of R.
766: // *
767: this._dormbr.Run("P", "L", "N", N, NRHS, N
768: , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref B, offset_b, LDB, ref WORK, NWORK + o_work
769: , LWORK - NWORK + 1, ref INFO);
770: // *
771: }
772: else
773: {
774: if (N >= MNTHR && LWORK >= 4 * M + M * M + Math.Max(M, Math.Max(2 * M - 4, Math.Max(NRHS, Math.Max(N - 3 * M, WLALSD)))))
775: {
776: // *
777: // * Path 2a - underdetermined, with many more columns than rows
778: // * and sufficient workspace for an efficient algorithm.
779: // *
780: LDWORK = M;
781: if (LWORK >= Math.Max(4 * M + M * LDA + Math.Max(M, Math.Max(2 * M - 4, Math.Max(NRHS, N - 3 * M))), Math.Max(M * LDA + M + M * NRHS, 4 * M + M * LDA + WLALSD))) LDWORK = LDA;
782: ITAU = 1;
783: NWORK = M + 1;
784: // *
785: // * Compute A=L*Q.
786: // * (Workspace: need 2*M, prefer M+M*NB)
787: // *
788: this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, NWORK + o_work
789: , LWORK - NWORK + 1, ref INFO);
790: IL = NWORK;
791: // *
792: // * Copy L to WORK(IL), zeroing out above its diagonal.
793: // *
794: this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref WORK, IL + o_work
795: , LDWORK);
796: this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref WORK, IL + LDWORK + o_work
797: , LDWORK);
798: IE = IL + LDWORK * M;
799: ITAUQ = IE + M;
800: ITAUP = ITAUQ + M;
801: NWORK = ITAUP + M;
802: // *
803: // * Bidiagonalize L in WORK(IL).
804: // * (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB)
805: // *
806: this._dgebrd.Run(M, M, ref WORK, IL + o_work, LDWORK, ref S, offset_s, ref WORK, IE + o_work
807: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref INFO);
808: // *
809: // * Multiply B by transpose of left bidiagonalizing vectors of L.
810: // * (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)
811: // *
812: this._dormbr.Run("Q", "L", "T", M, NRHS, M
813: , ref WORK, IL + o_work, LDWORK, WORK, ITAUQ + o_work, ref B, offset_b, LDB, ref WORK, NWORK + o_work
814: , LWORK - NWORK + 1, ref INFO);
815: // *
816: // * Solve the bidiagonal least squares problem.
817: // *
818: this._dlalsd.Run("U", SMLSIZ, M, NRHS, ref S, offset_s, ref WORK, IE + o_work
819: , ref B, offset_b, LDB, RCOND, ref RANK, ref WORK, NWORK + o_work, ref IWORK, offset_iwork
820: , ref INFO);
821: if (INFO != 0)
822: {
823: goto LABEL10;
824: }
825: // *
826: // * Multiply B by right bidiagonalizing vectors of L.
827: // *
828: this._dormbr.Run("P", "L", "N", M, NRHS, M
829: , ref WORK, IL + o_work, LDWORK, WORK, ITAUP + o_work, ref B, offset_b, LDB, ref WORK, NWORK + o_work
830: , LWORK - NWORK + 1, ref INFO);
831: // *
832: // * Zero out below first M rows of B.
833: // *
834: this._dlaset.Run("F", N - M, NRHS, ZERO, ZERO, ref B, M + 1+1 * LDB + o_b
835: , LDB);
836: NWORK = ITAU + M;
837: // *
838: // * Multiply transpose(Q) by B.
839: // * (Workspace: need M+NRHS, prefer M+NRHS*NB)
840: // *
841: this._dormlq.Run("L", "T", N, NRHS, M, ref A, offset_a
842: , LDA, WORK, ITAU + o_work, ref B, offset_b, LDB, ref WORK, NWORK + o_work, LWORK - NWORK + 1
843: , ref INFO);
844: // *
845: }
846: else
847: {
848: // *
849: // * Path 2 - remaining underdetermined cases.
850: // *
851: IE = 1;
852: ITAUQ = IE + M;
853: ITAUP = ITAUQ + M;
854: NWORK = ITAUP + M;
855: // *
856: // * Bidiagonalize A.
857: // * (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
858: // *
859: this._dgebrd.Run(M, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
860: , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref INFO);
861: // *
862: // * Multiply B by transpose of left bidiagonalizing vectors.
863: // * (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB)
864: // *
865: this._dormbr.Run("Q", "L", "T", M, NRHS, N
866: , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref B, offset_b, LDB, ref WORK, NWORK + o_work
867: , LWORK - NWORK + 1, ref INFO);
868: // *
869: // * Solve the bidiagonal least squares problem.
870: // *
871: this._dlalsd.Run("L", SMLSIZ, M, NRHS, ref S, offset_s, ref WORK, IE + o_work
872: , ref B, offset_b, LDB, RCOND, ref RANK, ref WORK, NWORK + o_work, ref IWORK, offset_iwork
873: , ref INFO);
874: if (INFO != 0)
875: {
876: goto LABEL10;
877: }
878: // *
879: // * Multiply B by right bidiagonalizing vectors of A.
880: // *
881: this._dormbr.Run("P", "L", "N", N, NRHS, M
882: , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref B, offset_b, LDB, ref WORK, NWORK + o_work
883: , LWORK - NWORK + 1, ref INFO);
884: // *
885: }
886: }
887: // *
888: // * Undo scaling.
889: // *
890: if (IASCL == 1)
891: {
892: this._dlascl.Run("G", 0, 0, ANRM, SMLNUM, N
893: , NRHS, ref B, offset_b, LDB, ref INFO);
894: this._dlascl.Run("G", 0, 0, SMLNUM, ANRM, MINMN
895: , 1, ref S, offset_s, MINMN, ref INFO);
896: }
897: else
898: {
899: if (IASCL == 2)
900: {
901: this._dlascl.Run("G", 0, 0, ANRM, BIGNUM, N
902: , NRHS, ref B, offset_b, LDB, ref INFO);
903: this._dlascl.Run("G", 0, 0, BIGNUM, ANRM, MINMN
904: , 1, ref S, offset_s, MINMN, ref INFO);
905: }
906: }
907: if (IBSCL == 1)
908: {
909: this._dlascl.Run("G", 0, 0, SMLNUM, BNRM, N
910: , NRHS, ref B, offset_b, LDB, ref INFO);
911: }
912: else
913: {
914: if (IBSCL == 2)
915: {
916: this._dlascl.Run("G", 0, 0, BIGNUM, BNRM, N
917: , NRHS, ref B, offset_b, LDB, ref INFO);
918: }
919: }
920: // *
921: LABEL10:;
922: WORK[1 + o_work] = MAXWRK;
923: return;
924: // *
925: // * End of DGELSD
926: // *
927:
928: #endregion
929:
930: }
931: }
932: }