1:  #region Translated by Jose Antonio De Santiago-Castillo.
2:
3:  //Translated by Jose Antonio De Santiago-Castillo.
4:  //E-mail:JAntonioDeSantiago@gmail.com
5:  //Web: www.DotNumerics.com
6:  //
7:  //Fortran to C# Translation.
8:  //Translated by:
9:  //F2CSharp Version 0.71 (November 10, 2009)
10:  //Code Optimizations: None
11:  //
12:  #endregion
13:
14:  using System;
15:  using DotNumerics.FortranLibrary;
16:
17:  namespace DotNumerics.CSLapack
18:  {
19:      /// <summary>
20:      /// -- LAPACK routine (version 3.1) --
21:      /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
22:      /// November 2006
23:      /// Purpose
24:      /// =======
25:      ///
26:      /// DORGBR generates one of the real orthogonal matrices Q or P**T
27:      /// determined by DGEBRD when reducing a real matrix A to bidiagonal
28:      /// form: A = Q * B * P**T.  Q and P**T are defined as products of
29:      /// elementary reflectors H(i) or G(i) respectively.
30:      ///
31:      /// If VECT = 'Q', A is assumed to have been an M-by-K matrix, and Q
32:      /// is of order M:
33:      /// if m .GE. k, Q = H(1) H(2) . . . H(k) and DORGBR returns the first n
34:      /// columns of Q, where m .GE. n .GE. k;
35:      /// if m .LT. k, Q = H(1) H(2) . . . H(m-1) and DORGBR returns Q as an
36:      /// M-by-M matrix.
37:      ///
38:      /// If VECT = 'P', A is assumed to have been a K-by-N matrix, and P**T
39:      /// is of order N:
40:      /// if k .LT. n, P**T = G(k) . . . G(2) G(1) and DORGBR returns the first m
41:      /// rows of P**T, where n .GE. m .GE. k;
42:      /// if k .GE. n, P**T = G(n-1) . . . G(2) G(1) and DORGBR returns P**T as
43:      /// an N-by-N matrix.
44:      ///
45:      ///</summary>
46:      public class DORGBR
47:      {
48:
49:
50:          #region Dependencies
51:
52:          LSAME _lsame; ILAENV _ilaenv; DORGLQ _dorglq; DORGQR _dorgqr; XERBLA _xerbla;
53:
54:          #endregion
55:
56:
57:          #region Fields
58:
59:          const double ZERO = 0.0E+0; const double ONE = 1.0E+0; bool LQUERY = false; bool WANTQ = false; int I = 0; int IINFO = 0;
60:          int J = 0;int LWKOPT = 0; int MN = 0; int NB = 0;
61:
62:          #endregion
63:
64:          public DORGBR(LSAME lsame, ILAENV ilaenv, DORGLQ dorglq, DORGQR dorgqr, XERBLA xerbla)
65:          {
66:
67:
68:              #region Set Dependencies
69:
70:              this._lsame = lsame; this._ilaenv = ilaenv; this._dorglq = dorglq; this._dorgqr = dorgqr; this._xerbla = xerbla;
71:
72:              #endregion
73:
74:          }
75:
76:          public DORGBR()
77:          {
78:
79:
80:              #region Dependencies (Initialization)
81:
82:              LSAME lsame = new LSAME();
83:              IEEECK ieeeck = new IEEECK();
84:              IPARMQ iparmq = new IPARMQ();
85:              DCOPY dcopy = new DCOPY();
86:              XERBLA xerbla = new XERBLA();
87:              DSCAL dscal = new DSCAL();
88:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
89:              DGEMM dgemm = new DGEMM(lsame, xerbla);
90:              DTRMM dtrmm = new DTRMM(lsame, xerbla);
91:              DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
92:              DGEMV dgemv = new DGEMV(lsame, xerbla);
93:              DTRMV dtrmv = new DTRMV(lsame, xerbla);
94:              DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);
95:              DGER dger = new DGER(xerbla);
96:              DLARF dlarf = new DLARF(dgemv, dger, lsame);
97:              DORGL2 dorgl2 = new DORGL2(dlarf, dscal, xerbla);
98:              DORGLQ dorglq = new DORGLQ(dlarfb, dlarft, dorgl2, xerbla, ilaenv);
99:              DORG2R dorg2r = new DORG2R(dlarf, dscal, xerbla);
100:              DORGQR dorgqr = new DORGQR(dlarfb, dlarft, dorg2r, xerbla, ilaenv);
101:
102:              #endregion
103:
104:
105:              #region Set Dependencies
106:
107:              this._lsame = lsame; this._ilaenv = ilaenv; this._dorglq = dorglq; this._dorgqr = dorgqr; this._xerbla = xerbla;
108:
109:              #endregion
110:
111:          }
112:          /// <summary>
113:          /// Purpose
114:          /// =======
115:          ///
116:          /// DORGBR generates one of the real orthogonal matrices Q or P**T
117:          /// determined by DGEBRD when reducing a real matrix A to bidiagonal
118:          /// form: A = Q * B * P**T.  Q and P**T are defined as products of
119:          /// elementary reflectors H(i) or G(i) respectively.
120:          ///
121:          /// If VECT = 'Q', A is assumed to have been an M-by-K matrix, and Q
122:          /// is of order M:
123:          /// if m .GE. k, Q = H(1) H(2) . . . H(k) and DORGBR returns the first n
124:          /// columns of Q, where m .GE. n .GE. k;
125:          /// if m .LT. k, Q = H(1) H(2) . . . H(m-1) and DORGBR returns Q as an
126:          /// M-by-M matrix.
127:          ///
128:          /// If VECT = 'P', A is assumed to have been a K-by-N matrix, and P**T
129:          /// is of order N:
130:          /// if k .LT. n, P**T = G(k) . . . G(2) G(1) and DORGBR returns the first m
131:          /// rows of P**T, where n .GE. m .GE. k;
132:          /// if k .GE. n, P**T = G(n-1) . . . G(2) G(1) and DORGBR returns P**T as
133:          /// an N-by-N matrix.
134:          ///
135:          ///</summary>
136:          /// <param name="VECT">
137:          /// (input) CHARACTER*1
138:          /// Specifies whether the matrix Q or the matrix P**T is
139:          /// required, as defined in the transformation applied by DGEBRD:
140:          /// = 'Q':  generate Q;
141:          /// = 'P':  generate P**T.
142:          ///</param>
143:          /// <param name="M">
144:          /// (input) INTEGER
145:          /// The number of rows of the matrix Q or P**T to be returned.
146:          /// M .GE. 0.
147:          ///</param>
148:          /// <param name="N">
149:          /// (input) INTEGER
150:          /// The number of columns of the matrix Q or P**T to be returned.
151:          /// N .GE. 0.
152:          /// If VECT = 'Q', M .GE. N .GE. min(M,K);
153:          /// if VECT = 'P', N .GE. M .GE. min(N,K).
154:          ///</param>
155:          /// <param name="K">
156:          /// (input) INTEGER
157:          /// If VECT = 'Q', the number of columns in the original M-by-K
158:          /// matrix reduced by DGEBRD.
159:          /// If VECT = 'P', the number of rows in the original K-by-N
160:          /// matrix reduced by DGEBRD.
161:          /// K .GE. 0.
162:          ///</param>
163:          /// <param name="A">
164:          /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
165:          /// On entry, the vectors which define the elementary reflectors,
166:          /// as returned by DGEBRD.
167:          /// On exit, the M-by-N matrix Q or P**T.
168:          ///</param>
169:          /// <param name="LDA">
170:          /// (input) INTEGER
171:          /// The leading dimension of the array A. LDA .GE. max(1,M).
172:          ///</param>
173:          /// <param name="TAU">
174:          /// (input) DOUBLE PRECISION array, dimension
175:          /// (min(M,K)) if VECT = 'Q'
176:          /// (min(N,K)) if VECT = 'P'
177:          /// TAU(i) must contain the scalar factor of the elementary
178:          /// reflector H(i) or G(i), which determines Q or P**T, as
179:          /// returned by DGEBRD in its array argument TAUQ or TAUP.
180:          ///</param>
181:          /// <param name="WORK">
182:          /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
183:          /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
184:          ///</param>
185:          /// <param name="LWORK">
186:          /// (input) INTEGER
187:          /// The dimension of the array WORK. LWORK .GE. max(1,min(M,N)).
188:          /// For optimum performance LWORK .GE. min(M,N)*NB, where NB
189:          /// is the optimal blocksize.
190:          ///
191:          /// If LWORK = -1, then a workspace query is assumed; the routine
192:          /// only calculates the optimal size of the WORK array, returns
193:          /// this value as the first entry of the WORK array, and no error
194:          /// message related to LWORK is issued by XERBLA.
195:          ///</param>
196:          /// <param name="INFO">
197:          /// (output) INTEGER
198:          /// = 0:  successful exit
199:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value
200:          ///</param>
201:          public void Run(string VECT, int M, int N, int K, ref double[] A, int offset_a, int LDA
202:                           , double[] TAU, int offset_tau, ref double[] WORK, int offset_work, int LWORK, ref int INFO)
203:          {
204:
205:              #region Array Index Correction
206:
207:               int o_a = -1 - LDA + offset_a;  int o_tau = -1 + offset_tau;  int o_work = -1 + offset_work;
208:
209:              #endregion
210:
211:
212:              #region Strings
213:
214:              VECT = VECT.Substring(0, 1);
215:
216:              #endregion
217:
218:
219:              #region Prolog
220:
221:              // *
222:              // *  -- LAPACK routine (version 3.1) --
223:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
224:              // *     November 2006
225:              // *
226:              // *     .. Scalar Arguments ..
227:              // *     ..
228:              // *     .. Array Arguments ..
229:              // *     ..
230:              // *
231:              // *  Purpose
232:              // *  =======
233:              // *
234:              // *  DORGBR generates one of the real orthogonal matrices Q or P**T
235:              // *  determined by DGEBRD when reducing a real matrix A to bidiagonal
236:              // *  form: A = Q * B * P**T.  Q and P**T are defined as products of
237:              // *  elementary reflectors H(i) or G(i) respectively.
238:              // *
239:              // *  If VECT = 'Q', A is assumed to have been an M-by-K matrix, and Q
240:              // *  is of order M:
241:              // *  if m >= k, Q = H(1) H(2) . . . H(k) and DORGBR returns the first n
242:              // *  columns of Q, where m >= n >= k;
243:              // *  if m < k, Q = H(1) H(2) . . . H(m-1) and DORGBR returns Q as an
244:              // *  M-by-M matrix.
245:              // *
246:              // *  If VECT = 'P', A is assumed to have been a K-by-N matrix, and P**T
247:              // *  is of order N:
248:              // *  if k < n, P**T = G(k) . . . G(2) G(1) and DORGBR returns the first m
249:              // *  rows of P**T, where n >= m >= k;
250:              // *  if k >= n, P**T = G(n-1) . . . G(2) G(1) and DORGBR returns P**T as
251:              // *  an N-by-N matrix.
252:              // *
253:              // *  Arguments
254:              // *  =========
255:              // *
256:              // *  VECT    (input) CHARACTER*1
257:              // *          Specifies whether the matrix Q or the matrix P**T is
258:              // *          required, as defined in the transformation applied by DGEBRD:
259:              // *          = 'Q':  generate Q;
260:              // *          = 'P':  generate P**T.
261:              // *
262:              // *  M       (input) INTEGER
263:              // *          The number of rows of the matrix Q or P**T to be returned.
264:              // *          M >= 0.
265:              // *
266:              // *  N       (input) INTEGER
267:              // *          The number of columns of the matrix Q or P**T to be returned.
268:              // *          N >= 0.
269:              // *          If VECT = 'Q', M >= N >= min(M,K);
270:              // *          if VECT = 'P', N >= M >= min(N,K).
271:              // *
272:              // *  K       (input) INTEGER
273:              // *          If VECT = 'Q', the number of columns in the original M-by-K
274:              // *          matrix reduced by DGEBRD.
275:              // *          If VECT = 'P', the number of rows in the original K-by-N
276:              // *          matrix reduced by DGEBRD.
277:              // *          K >= 0.
278:              // *
279:              // *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
280:              // *          On entry, the vectors which define the elementary reflectors,
281:              // *          as returned by DGEBRD.
282:              // *          On exit, the M-by-N matrix Q or P**T.
283:              // *
284:              // *  LDA     (input) INTEGER
285:              // *          The leading dimension of the array A. LDA >= max(1,M).
286:              // *
287:              // *  TAU     (input) DOUBLE PRECISION array, dimension
288:              // *                                (min(M,K)) if VECT = 'Q'
289:              // *                                (min(N,K)) if VECT = 'P'
290:              // *          TAU(i) must contain the scalar factor of the elementary
291:              // *          reflector H(i) or G(i), which determines Q or P**T, as
292:              // *          returned by DGEBRD in its array argument TAUQ or TAUP.
293:              // *
294:              // *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
295:              // *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
296:              // *
297:              // *  LWORK   (input) INTEGER
298:              // *          The dimension of the array WORK. LWORK >= max(1,min(M,N)).
299:              // *          For optimum performance LWORK >= min(M,N)*NB, where NB
300:              // *          is the optimal blocksize.
301:              // *
302:              // *          If LWORK = -1, then a workspace query is assumed; the routine
303:              // *          only calculates the optimal size of the WORK array, returns
304:              // *          this value as the first entry of the WORK array, and no error
305:              // *          message related to LWORK is issued by XERBLA.
306:              // *
307:              // *  INFO    (output) INTEGER
308:              // *          = 0:  successful exit
309:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value
310:              // *
311:              // *  =====================================================================
312:              // *
313:              // *     .. Parameters ..
314:              // *     ..
315:              // *     .. Local Scalars ..
316:              // *     ..
317:              // *     .. External Functions ..
318:              // *     ..
319:              // *     .. External Subroutines ..
320:              // *     ..
321:              // *     .. Intrinsic Functions ..
322:              //      INTRINSIC          MAX, MIN;
323:              // *     ..
324:              // *     .. Executable Statements ..
325:              // *
326:              // *     Test the input arguments
327:              // *
328:
329:              #endregion
330:
331:
332:              #region Body
333:
334:              INFO = 0;
335:              WANTQ = this._lsame.Run(VECT, "Q");
336:              MN = Math.Min(M, N);
337:              LQUERY = (LWORK ==  - 1);
338:              if (!WANTQ && !this._lsame.Run(VECT, "P"))
339:              {
340:                  INFO =  - 1;
341:              }
342:              else
343:              {
344:                  if (M < 0)
345:                  {
346:                      INFO =  - 2;
347:                  }
348:                  else
349:                  {
350:                      if (N < 0 || (WANTQ && (N > M || N < Math.Min(M, K))) || (!WANTQ && (M > N || M < Math.Min(N, K))))
351:                      {
352:                          INFO =  - 3;
353:                      }
354:                      else
355:                      {
356:                          if (K < 0)
357:                          {
358:                              INFO =  - 4;
359:                          }
360:                          else
361:                          {
362:                              if (LDA < Math.Max(1, M))
363:                              {
364:                                  INFO =  - 6;
365:                              }
366:                              else
367:                              {
368:                                  if (LWORK < Math.Max(1, MN) && !LQUERY)
369:                                  {
370:                                      INFO =  - 9;
371:                                  }
372:                              }
373:                          }
374:                      }
375:                  }
376:              }
377:              // *
378:              if (INFO == 0)
379:              {
380:                  if (WANTQ)
381:                  {
382:                      NB = this._ilaenv.Run(1, "DORGQR", " ", M, N, K,  - 1);
383:                  }
384:                  else
385:                  {
386:                      NB = this._ilaenv.Run(1, "DORGLQ", " ", M, N, K,  - 1);
387:                  }
388:                  LWKOPT = Math.Max(1, MN) * NB;
389:                  WORK[1 + o_work] = LWKOPT;
390:              }
391:              // *
392:              if (INFO != 0)
393:              {
394:                  this._xerbla.Run("DORGBR",  - INFO);
395:                  return;
396:              }
397:              else
398:              {
399:                  if (LQUERY)
400:                  {
401:                      return;
402:                  }
403:              }
404:              // *
405:              // *     Quick return if possible
406:              // *
407:              if (M == 0 || N == 0)
408:              {
409:                  WORK[1 + o_work] = 1;
410:                  return;
411:              }
412:              // *
413:              if (WANTQ)
414:              {
415:                  // *
416:                  // *        Form Q, determined by a call to DGEBRD to reduce an m-by-k
417:                  // *        matrix
418:                  // *
419:                  if (M >= K)
420:                  {
421:                      // *
422:                      // *           If m >= k, assume m >= n >= k
423:                      // *
424:                      this._dorgqr.Run(M, N, K, ref A, offset_a, LDA, TAU, offset_tau
425:                                       , ref WORK, offset_work, LWORK, ref IINFO);
426:                      // *
427:                  }
428:                  else
429:                  {
430:                      // *
431:                      // *           If m < k, assume m = n
432:                      // *
433:                      // *           Shift the vectors which define the elementary reflectors one
434:                      // *           column to the right, and set the first row and column of Q
435:                      // *           to those of the unit matrix
436:                      // *
437:                      for (J = M; J >= 2; J +=  - 1)
438:                      {
439:                          A[1+J * LDA + o_a] = ZERO;
440:                          for (I = J + 1; I <= M; I++)
441:                          {
442:                              A[I+J * LDA + o_a] = A[I+(J - 1) * LDA + o_a];
443:                          }
444:                      }
445:                      A[1+1 * LDA + o_a] = ONE;
446:                      for (I = 2; I <= M; I++)
447:                      {
448:                          A[I+1 * LDA + o_a] = ZERO;
449:                      }
450:                      if (M > 1)
451:                      {
452:                          // *
453:                          // *              Form Q(2:m,2:m)
454:                          // *
455:                          this._dorgqr.Run(M - 1, M - 1, M - 1, ref A, 2+2 * LDA + o_a, LDA, TAU, offset_tau
456:                                           , ref WORK, offset_work, LWORK, ref IINFO);
457:                      }
458:                  }
459:              }
460:              else
461:              {
462:                  // *
463:                  // *        Form P', determined by a call to DGEBRD to reduce a k-by-n
464:                  // *        matrix
465:                  // *
466:                  if (K < N)
467:                  {
468:                      // *
469:                      // *           If k < n, assume k <= m <= n
470:                      // *
471:                      this._dorglq.Run(M, N, K, ref A, offset_a, LDA, TAU, offset_tau
472:                                       , ref WORK, offset_work, LWORK, ref IINFO);
473:                      // *
474:                  }
475:                  else
476:                  {
477:                      // *
478:                      // *           If k >= n, assume m = n
479:                      // *
480:                      // *           Shift the vectors which define the elementary reflectors one
481:                      // *           row downward, and set the first row and column of P' to
482:                      // *           those of the unit matrix
483:                      // *
484:                      A[1+1 * LDA + o_a] = ONE;
485:                      for (I = 2; I <= N; I++)
486:                      {
487:                          A[I+1 * LDA + o_a] = ZERO;
488:                      }
489:                      for (J = 2; J <= N; J++)
490:                      {
491:                          for (I = J - 1; I >= 2; I +=  - 1)
492:                          {
493:                              A[I+J * LDA + o_a] = A[I - 1+J * LDA + o_a];
494:                          }
495:                          A[1+J * LDA + o_a] = ZERO;
496:                      }
497:                      if (N > 1)
498:                      {
499:                          // *
500:                          // *              Form P'(2:n,2:n)
501:                          // *
502:                          this._dorglq.Run(N - 1, N - 1, N - 1, ref A, 2+2 * LDA + o_a, LDA, TAU, offset_tau
503:                                           , ref WORK, offset_work, LWORK, ref IINFO);
504:                      }
505:                  }
506:              }
507:              WORK[1 + o_work] = LWKOPT;
508:              return;
509:              // *
510:              // *     End of DORGBR
511:              // *
512:
513:              #endregion
514:
515:          }
516:      }
517:  }