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: /// DGBTRF computes an LU factorization of a real m-by-n band matrix A
27: /// using partial pivoting with row interchanges.
28: ///
29: /// This is the blocked version of the algorithm, calling Level 3 BLAS.
30: ///
31: ///</summary>
32: public class DGBTRF
33: {
34:
35:
36: #region Dependencies
37:
38: IDAMAX _idamax; ILAENV _ilaenv; DCOPY _dcopy; DGBTF2 _dgbtf2; DGEMM _dgemm; DGER _dger; DLASWP _dlaswp; DSCAL _dscal;
39: DSWAP _dswap;DTRSM _dtrsm; XERBLA _xerbla;
40:
41: #endregion
42:
43:
44: #region Fields
45:
46: const double ONE = 1.0E+0; const double ZERO = 0.0E+0; const int NBMAX = 64; const int LDWORK = NBMAX + 1; int I = 0;
47: int I2 = 0;int I3 = 0; int II = 0; int IP = 0; int J = 0; int J2 = 0; int J3 = 0; int JB = 0; int JJ = 0; int JM = 0;
48: int JP = 0;int JU = 0; int K2 = 0; int KM = 0; int KV = 0; int NB = 0; int NW = 0; double TEMP = 0;
49: double[] WORK13 = new double[LDWORK * NBMAX]; int offset_work13 = 0; int o_work13 = -1 - LDWORK;
50: double[] WORK31 = new double[LDWORK * NBMAX]; int offset_work31 = 0; int o_work31 = -1 - LDWORK;
51:
52: #endregion
53:
54: public DGBTRF(IDAMAX idamax, ILAENV ilaenv, DCOPY dcopy, DGBTF2 dgbtf2, DGEMM dgemm, DGER dger, DLASWP dlaswp, DSCAL dscal, DSWAP dswap, DTRSM dtrsm
55: , XERBLA xerbla)
56: {
57:
58:
59: #region Set Dependencies
60:
61: this._idamax = idamax; this._ilaenv = ilaenv; this._dcopy = dcopy; this._dgbtf2 = dgbtf2; this._dgemm = dgemm;
62: this._dger = dger;this._dlaswp = dlaswp; this._dscal = dscal; this._dswap = dswap; this._dtrsm = dtrsm;
63: this._xerbla = xerbla;
64:
65: #endregion
66:
67: }
68:
69: public DGBTRF()
70: {
71:
72:
73: #region Dependencies (Initialization)
74:
75: IDAMAX idamax = new IDAMAX();
76: IEEECK ieeeck = new IEEECK();
77: IPARMQ iparmq = new IPARMQ();
78: DCOPY dcopy = new DCOPY();
79: XERBLA xerbla = new XERBLA();
80: DSCAL dscal = new DSCAL();
81: DSWAP dswap = new DSWAP();
82: LSAME lsame = new LSAME();
83: DLASWP dlaswp = new DLASWP();
84: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
85: DGER dger = new DGER(xerbla);
86: DGBTF2 dgbtf2 = new DGBTF2(idamax, dger, dscal, dswap, xerbla);
87: DGEMM dgemm = new DGEMM(lsame, xerbla);
88: DTRSM dtrsm = new DTRSM(lsame, xerbla);
89:
90: #endregion
91:
92:
93: #region Set Dependencies
94:
95: this._idamax = idamax; this._ilaenv = ilaenv; this._dcopy = dcopy; this._dgbtf2 = dgbtf2; this._dgemm = dgemm;
96: this._dger = dger;this._dlaswp = dlaswp; this._dscal = dscal; this._dswap = dswap; this._dtrsm = dtrsm;
97: this._xerbla = xerbla;
98:
99: #endregion
100:
101: }
102: /// <summary>
103: /// Purpose
104: /// =======
105: ///
106: /// DGBTRF computes an LU factorization of a real m-by-n band matrix A
107: /// using partial pivoting with row interchanges.
108: ///
109: /// This is the blocked version of the algorithm, calling Level 3 BLAS.
110: ///
111: ///</summary>
112: /// <param name="M">
113: /// (input) INTEGER
114: /// The number of rows of the matrix A. M .GE. 0.
115: ///</param>
116: /// <param name="N">
117: /// (input) INTEGER
118: /// The number of columns of the matrix A. N .GE. 0.
119: ///</param>
120: /// <param name="KL">
121: /// (input) INTEGER
122: /// The number of subdiagonals within the band of A. KL .GE. 0.
123: ///</param>
124: /// <param name="KU">
125: /// (input) INTEGER
126: /// The number of superdiagonals within the band of A. KU .GE. 0.
127: ///</param>
128: /// <param name="AB">
129: /// (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
130: /// On entry, the matrix A in band storage, in rows KL+1 to
131: /// 2*KL+KU+1; rows 1 to KL of the array need not be set.
132: /// The j-th column of A is stored in the j-th column of the
133: /// array AB as follows:
134: /// AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku).LE.i.LE.min(m,j+kl)
135: ///
136: /// On exit, details of the factorization: U is stored as an
137: /// upper triangular band matrix with KL+KU superdiagonals in
138: /// rows 1 to KL+KU+1, and the multipliers used during the
139: /// factorization are stored in rows KL+KU+2 to 2*KL+KU+1.
140: /// See below for further details.
141: ///</param>
142: /// <param name="LDAB">
143: /// (input) INTEGER
144: /// The leading dimension of the array AB. LDAB .GE. 2*KL+KU+1.
145: ///</param>
146: /// <param name="IPIV">
147: /// (output) INTEGER array, dimension (min(M,N))
148: /// The pivot indices; for 1 .LE. i .LE. min(M,N), row i of the
149: /// matrix was interchanged with row IPIV(i).
150: ///</param>
151: /// <param name="INFO">
152: /// (output) INTEGER
153: /// = 0: successful exit
154: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value
155: /// .GT. 0: if INFO = +i, U(i,i) is exactly zero. The factorization
156: /// has been completed, but the factor U is exactly
157: /// singular, and division by zero will occur if it is used
158: /// to solve a system of equations.
159: ///</param>
160: public void Run(int M, int N, int KL, int KU, ref double[] AB, int offset_ab, int LDAB
161: , ref int[] IPIV, int offset_ipiv, ref int INFO)
162: {
163:
164: #region Array Index Correction
165:
166: int o_ab = -1 - LDAB + offset_ab; int o_ipiv = -1 + offset_ipiv;
167:
168: #endregion
169:
170:
171: #region Prolog
172:
173: // *
174: // * -- LAPACK routine (version 3.1) --
175: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
176: // * November 2006
177: // *
178: // * .. Scalar Arguments ..
179: // * ..
180: // * .. Array Arguments ..
181: // * ..
182: // *
183: // * Purpose
184: // * =======
185: // *
186: // * DGBTRF computes an LU factorization of a real m-by-n band matrix A
187: // * using partial pivoting with row interchanges.
188: // *
189: // * This is the blocked version of the algorithm, calling Level 3 BLAS.
190: // *
191: // * Arguments
192: // * =========
193: // *
194: // * M (input) INTEGER
195: // * The number of rows of the matrix A. M >= 0.
196: // *
197: // * N (input) INTEGER
198: // * The number of columns of the matrix A. N >= 0.
199: // *
200: // * KL (input) INTEGER
201: // * The number of subdiagonals within the band of A. KL >= 0.
202: // *
203: // * KU (input) INTEGER
204: // * The number of superdiagonals within the band of A. KU >= 0.
205: // *
206: // * AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
207: // * On entry, the matrix A in band storage, in rows KL+1 to
208: // * 2*KL+KU+1; rows 1 to KL of the array need not be set.
209: // * The j-th column of A is stored in the j-th column of the
210: // * array AB as follows:
211: // * AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl)
212: // *
213: // * On exit, details of the factorization: U is stored as an
214: // * upper triangular band matrix with KL+KU superdiagonals in
215: // * rows 1 to KL+KU+1, and the multipliers used during the
216: // * factorization are stored in rows KL+KU+2 to 2*KL+KU+1.
217: // * See below for further details.
218: // *
219: // * LDAB (input) INTEGER
220: // * The leading dimension of the array AB. LDAB >= 2*KL+KU+1.
221: // *
222: // * IPIV (output) INTEGER array, dimension (min(M,N))
223: // * The pivot indices; for 1 <= i <= min(M,N), row i of the
224: // * matrix was interchanged with row IPIV(i).
225: // *
226: // * INFO (output) INTEGER
227: // * = 0: successful exit
228: // * < 0: if INFO = -i, the i-th argument had an illegal value
229: // * > 0: if INFO = +i, U(i,i) is exactly zero. The factorization
230: // * has been completed, but the factor U is exactly
231: // * singular, and division by zero will occur if it is used
232: // * to solve a system of equations.
233: // *
234: // * Further Details
235: // * ===============
236: // *
237: // * The band storage scheme is illustrated by the following example, when
238: // * M = N = 6, KL = 2, KU = 1:
239: // *
240: // * On entry: On exit:
241: // *
242: // * * * * + + + * * * u14 u25 u36
243: // * * * + + + + * * u13 u24 u35 u46
244: // * * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56
245: // * a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66
246: // * a21 a32 a43 a54 a65 * m21 m32 m43 m54 m65 *
247: // * a31 a42 a53 a64 * * m31 m42 m53 m64 * *
248: // *
249: // * Array elements marked * are not used by the routine; elements marked
250: // * + need not be set on entry, but are required by the routine to store
251: // * elements of U because of fill-in resulting from the row interchanges.
252: // *
253: // * =====================================================================
254: // *
255: // * .. Parameters ..
256: // * ..
257: // * .. Local Scalars ..
258: // * ..
259: // * .. Local Arrays ..
260: // * ..
261: // * .. External Functions ..
262: // * ..
263: // * .. External Subroutines ..
264: // * ..
265: // * .. Intrinsic Functions ..
266: // INTRINSIC MAX, MIN;
267: // * ..
268: // * .. Executable Statements ..
269: // *
270: // * KV is the number of superdiagonals in the factor U, allowing for
271: // * fill-in
272: // *
273:
274: #endregion
275:
276:
277: #region Body
278:
279: KV = KU + KL;
280: // *
281: // * Test the input parameters.
282: // *
283: INFO = 0;
284: if (M < 0)
285: {
286: INFO = - 1;
287: }
288: else
289: {
290: if (N < 0)
291: {
292: INFO = - 2;
293: }
294: else
295: {
296: if (KL < 0)
297: {
298: INFO = - 3;
299: }
300: else
301: {
302: if (KU < 0)
303: {
304: INFO = - 4;
305: }
306: else
307: {
308: if (LDAB < KL + KV + 1)
309: {
310: INFO = - 6;
311: }
312: }
313: }
314: }
315: }
316: if (INFO != 0)
317: {
318: this._xerbla.Run("DGBTRF", - INFO);
319: return;
320: }
321: // *
322: // * Quick return if possible
323: // *
324: if (M == 0 || N == 0) return;
325: // *
326: // * Determine the block size for this environment
327: // *
328: NB = this._ilaenv.Run(1, "DGBTRF", " ", M, N, KL, KU);
329: // *
330: // * The block size must not exceed the limit set by the size of the
331: // * local arrays WORK13 and WORK31.
332: // *
333: NB = Math.Min(NB, NBMAX);
334: // *
335: if (NB <= 1 || NB > KL)
336: {
337: // *
338: // * Use unblocked code
339: // *
340: this._dgbtf2.Run(M, N, KL, KU, ref AB, offset_ab, LDAB
341: , ref IPIV, offset_ipiv, ref INFO);
342: }
343: else
344: {
345: // *
346: // * Use blocked code
347: // *
348: // * Zero the superdiagonal elements of the work array WORK13
349: // *
350: for (J = 1; J <= NB; J++)
351: {
352: for (I = 1; I <= J - 1; I++)
353: {
354: WORK13[I+J * LDWORK + o_work13] = ZERO;
355: }
356: }
357: // *
358: // * Zero the subdiagonal elements of the work array WORK31
359: // *
360: for (J = 1; J <= NB; J++)
361: {
362: for (I = J + 1; I <= NB; I++)
363: {
364: WORK31[I+J * LDWORK + o_work31] = ZERO;
365: }
366: }
367: // *
368: // * Gaussian elimination with partial pivoting
369: // *
370: // * Set fill-in elements in columns KU+2 to KV to zero
371: // *
372: for (J = KU + 2; J <= Math.Min(KV, N); J++)
373: {
374: for (I = KV - J + 2; I <= KL; I++)
375: {
376: AB[I+J * LDAB + o_ab] = ZERO;
377: }
378: }
379: // *
380: // * JU is the index of the last column affected by the current
381: // * stage of the factorization
382: // *
383: JU = 1;
384: // *
385: for (J = 1; (NB >= 0) ? (J <= Math.Min(M, N)) : (J >= Math.Min(M, N)); J += NB)
386: {
387: JB = Math.Min(NB, Math.Min(M, N) - J + 1);
388: // *
389: // * The active part of the matrix is partitioned
390: // *
391: // * A11 A12 A13
392: // * A21 A22 A23
393: // * A31 A32 A33
394: // *
395: // * Here A11, A21 and A31 denote the current block of JB columns
396: // * which is about to be factorized. The number of rows in the
397: // * partitioning are JB, I2, I3 respectively, and the numbers
398: // * of columns are JB, J2, J3. The superdiagonal elements of A13
399: // * and the subdiagonal elements of A31 lie outside the band.
400: // *
401: I2 = Math.Min(KL - JB, M - J - JB + 1);
402: I3 = Math.Min(JB, M - J - KL + 1);
403: // *
404: // * J2 and J3 are computed after JU has been updated.
405: // *
406: // * Factorize the current block of JB columns
407: // *
408: for (JJ = J; JJ <= J + JB - 1; JJ++)
409: {
410: // *
411: // * Set fill-in elements in column JJ+KV to zero
412: // *
413: if (JJ + KV <= N)
414: {
415: for (I = 1; I <= KL; I++)
416: {
417: AB[I+(JJ + KV) * LDAB + o_ab] = ZERO;
418: }
419: }
420: // *
421: // * Find pivot and test for singularity. KM is the number of
422: // * subdiagonal elements in the current column.
423: // *
424: KM = Math.Min(KL, M - JJ);
425: JP = this._idamax.Run(KM + 1, AB, KV + 1+JJ * LDAB + o_ab, 1);
426: IPIV[JJ + o_ipiv] = JP + JJ - J;
427: if (AB[KV + JP+JJ * LDAB + o_ab] != ZERO)
428: {
429: JU = Math.Max(JU, Math.Min(JJ + KU + JP - 1, N));
430: if (JP != 1)
431: {
432: // *
433: // * Apply interchange to columns J to J+JB-1
434: // *
435: if (JP + JJ - 1 < J + KL)
436: {
437: // *
438: this._dswap.Run(JB, ref AB, KV + 1 + JJ - J+J * LDAB + o_ab, LDAB - 1, ref AB, KV + JP + JJ - J+J * LDAB + o_ab, LDAB - 1);
439: }
440: else
441: {
442: // *
443: // * The interchange affects columns J to JJ-1 of A31
444: // * which are stored in the work array WORK31
445: // *
446: this._dswap.Run(JJ - J, ref AB, KV + 1 + JJ - J+J * LDAB + o_ab, LDAB - 1, ref WORK31, JP + JJ - J - KL+1 * LDWORK + o_work31, LDWORK);
447: this._dswap.Run(J + JB - JJ, ref AB, KV + 1+JJ * LDAB + o_ab, LDAB - 1, ref AB, KV + JP+JJ * LDAB + o_ab, LDAB - 1);
448: }
449: }
450: // *
451: // * Compute multipliers
452: // *
453: this._dscal.Run(KM, ONE / AB[KV + 1+JJ * LDAB + o_ab], ref AB, KV + 2+JJ * LDAB + o_ab, 1);
454: // *
455: // * Update trailing submatrix within the band and within
456: // * the current block. JM is the index of the last column
457: // * which needs to be updated.
458: // *
459: JM = Math.Min(JU, J + JB - 1);
460: if (JM > JJ)
461: {
462: this._dger.Run(KM, JM - JJ, - ONE, AB, KV + 2+JJ * LDAB + o_ab, 1, AB, KV+(JJ + 1) * LDAB + o_ab
463: , LDAB - 1, ref AB, KV + 1+(JJ + 1) * LDAB + o_ab, LDAB - 1);
464: }
465: }
466: else
467: {
468: // *
469: // * If pivot is zero, set INFO to the index of the pivot
470: // * unless a zero pivot has already been found.
471: // *
472: if (INFO == 0) INFO = JJ;
473: }
474: // *
475: // * Copy current column of A31 into the work array WORK31
476: // *
477: NW = Math.Min(JJ - J + 1, I3);
478: if (NW > 0) this._dcopy.Run(NW, AB, KV + KL + 1 - JJ + J+JJ * LDAB + o_ab, 1, ref WORK31, 1+(JJ - J + 1) * LDWORK + o_work31, 1);
479: }
480: if (J + JB <= N)
481: {
482: // *
483: // * Apply the row interchanges to the other blocks.
484: // *
485: J2 = Math.Min(JU - J + 1, KV) - JB;
486: J3 = Math.Max(0, JU - J - KV + 1);
487: // *
488: // * Use DLASWP to apply the row interchanges to A12, A22, and
489: // * A32.
490: // *
491: this._dlaswp.Run(J2, ref AB, KV + 1 - JB+(J + JB) * LDAB + o_ab, LDAB - 1, 1, JB, IPIV, J + o_ipiv
492: , 1);
493: // *
494: // * Adjust the pivot indices.
495: // *
496: for (I = J; I <= J + JB - 1; I++)
497: {
498: IPIV[I + o_ipiv] = IPIV[I + o_ipiv] + J - 1;
499: }
500: // *
501: // * Apply the row interchanges to A13, A23, and A33
502: // * columnwise.
503: // *
504: K2 = J - 1 + JB + J2;
505: for (I = 1; I <= J3; I++)
506: {
507: JJ = K2 + I;
508: for (II = J + I - 1; II <= J + JB - 1; II++)
509: {
510: IP = IPIV[II + o_ipiv];
511: if (IP != II)
512: {
513: TEMP = AB[KV + 1 + II - JJ+JJ * LDAB + o_ab];
514: AB[KV + 1 + II - JJ+JJ * LDAB + o_ab] = AB[KV + 1 + IP - JJ+JJ * LDAB + o_ab];
515: AB[KV + 1 + IP - JJ+JJ * LDAB + o_ab] = TEMP;
516: }
517: }
518: }
519: // *
520: // * Update the relevant part of the trailing submatrix
521: // *
522: if (J2 > 0)
523: {
524: // *
525: // * Update A12
526: // *
527: this._dtrsm.Run("Left", "Lower", "No transpose", "Unit", JB, J2
528: , ONE, AB, KV + 1+J * LDAB + o_ab, LDAB - 1, ref AB, KV + 1 - JB+(J + JB) * LDAB + o_ab, LDAB - 1);
529: // *
530: if (I2 > 0)
531: {
532: // *
533: // * Update A22
534: // *
535: this._dgemm.Run("No transpose", "No transpose", I2, J2, JB, - ONE
536: , AB, KV + 1 + JB+J * LDAB + o_ab, LDAB - 1, AB, KV + 1 - JB+(J + JB) * LDAB + o_ab, LDAB - 1, ONE, ref AB, KV + 1+(J + JB) * LDAB + o_ab
537: , LDAB - 1);
538: }
539: // *
540: if (I3 > 0)
541: {
542: // *
543: // * Update A32
544: // *
545: this._dgemm.Run("No transpose", "No transpose", I3, J2, JB, - ONE
546: , WORK31, offset_work31, LDWORK, AB, KV + 1 - JB+(J + JB) * LDAB + o_ab, LDAB - 1, ONE, ref AB, KV + KL + 1 - JB+(J + JB) * LDAB + o_ab
547: , LDAB - 1);
548: }
549: }
550: // *
551: if (J3 > 0)
552: {
553: // *
554: // * Copy the lower triangle of A13 into the work array
555: // * WORK13
556: // *
557: for (JJ = 1; JJ <= J3; JJ++)
558: {
559: for (II = JJ; II <= JB; II++)
560: {
561: WORK13[II+JJ * LDWORK + o_work13] = AB[II - JJ + 1+(JJ + J + KV - 1) * LDAB + o_ab];
562: }
563: }
564: // *
565: // * Update A13 in the work array
566: // *
567: this._dtrsm.Run("Left", "Lower", "No transpose", "Unit", JB, J3
568: , ONE, AB, KV + 1+J * LDAB + o_ab, LDAB - 1, ref WORK13, offset_work13, LDWORK);
569: // *
570: if (I2 > 0)
571: {
572: // *
573: // * Update A23
574: // *
575: this._dgemm.Run("No transpose", "No transpose", I2, J3, JB, - ONE
576: , AB, KV + 1 + JB+J * LDAB + o_ab, LDAB - 1, WORK13, offset_work13, LDWORK, ONE, ref AB, 1 + JB+(J + KV) * LDAB + o_ab
577: , LDAB - 1);
578: }
579: // *
580: if (I3 > 0)
581: {
582: // *
583: // * Update A33
584: // *
585: this._dgemm.Run("No transpose", "No transpose", I3, J3, JB, - ONE
586: , WORK31, offset_work31, LDWORK, WORK13, offset_work13, LDWORK, ONE, ref AB, 1 + KL+(J + KV) * LDAB + o_ab
587: , LDAB - 1);
588: }
589: // *
590: // * Copy the lower triangle of A13 back into place
591: // *
592: for (JJ = 1; JJ <= J3; JJ++)
593: {
594: for (II = JJ; II <= JB; II++)
595: {
596: AB[II - JJ + 1+(JJ + J + KV - 1) * LDAB + o_ab] = WORK13[II+JJ * LDWORK + o_work13];
597: }
598: }
599: }
600: }
601: else
602: {
603: // *
604: // * Adjust the pivot indices.
605: // *
606: for (I = J; I <= J + JB - 1; I++)
607: {
608: IPIV[I + o_ipiv] = IPIV[I + o_ipiv] + J - 1;
609: }
610: }
611: // *
612: // * Partially undo the interchanges in the current block to
613: // * restore the upper triangular form of A31 and copy the upper
614: // * triangle of A31 back into place
615: // *
616: for (JJ = J + JB - 1; JJ >= J; JJ += - 1)
617: {
618: JP = IPIV[JJ + o_ipiv] - JJ + 1;
619: if (JP != 1)
620: {
621: // *
622: // * Apply interchange to columns J to JJ-1
623: // *
624: if (JP + JJ - 1 < J + KL)
625: {
626: // *
627: // * The interchange does not affect A31
628: // *
629: this._dswap.Run(JJ - J, ref AB, KV + 1 + JJ - J+J * LDAB + o_ab, LDAB - 1, ref AB, KV + JP + JJ - J+J * LDAB + o_ab, LDAB - 1);
630: }
631: else
632: {
633: // *
634: // * The interchange does affect A31
635: // *
636: this._dswap.Run(JJ - J, ref AB, KV + 1 + JJ - J+J * LDAB + o_ab, LDAB - 1, ref WORK31, JP + JJ - J - KL+1 * LDWORK + o_work31, LDWORK);
637: }
638: }
639: // *
640: // * Copy the current column of A31 back into place
641: // *
642: NW = Math.Min(I3, JJ - J + 1);
643: if (NW > 0) this._dcopy.Run(NW, WORK31, 1+(JJ - J + 1) * LDWORK + o_work31, 1, ref AB, KV + KL + 1 - JJ + J+JJ * LDAB + o_ab, 1);
644: }
645: }
646: }
647: // *
648: return;
649: // *
650: // * End of DGBTRF
651: // *
652:
653: #endregion
654:
655: }
656: }
657: }