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: /// DLAEDA computes the Z vector corresponding to the merge step in the
27: /// CURLVLth step of the merge process with TLVLS steps for the CURPBMth
28: /// problem.
29: ///
30: ///</summary>
31: public class DLAEDA
32: {
33:
34:
35: #region Dependencies
36:
37: DCOPY _dcopy; DGEMV _dgemv; DROT _drot; XERBLA _xerbla;
38:
39: #endregion
40:
41:
42: #region Fields
43:
44: const double ZERO = 0.0E0; const double HALF = 0.5E0; const double ONE = 1.0E0; int BSIZ1 = 0; int BSIZ2 = 0;
45: int CURR = 0;int I = 0; int K = 0; int MID = 0; int PSIZ1 = 0; int PSIZ2 = 0; int PTR = 0; int ZPTR1 = 0;
46:
47: #endregion
48:
49: public DLAEDA(DCOPY dcopy, DGEMV dgemv, DROT drot, XERBLA xerbla)
50: {
51:
52:
53: #region Set Dependencies
54:
55: this._dcopy = dcopy; this._dgemv = dgemv; this._drot = drot; this._xerbla = xerbla;
56:
57: #endregion
58:
59: }
60:
61: public DLAEDA()
62: {
63:
64:
65: #region Dependencies (Initialization)
66:
67: DCOPY dcopy = new DCOPY();
68: LSAME lsame = new LSAME();
69: XERBLA xerbla = new XERBLA();
70: DROT drot = new DROT();
71: DGEMV dgemv = new DGEMV(lsame, xerbla);
72:
73: #endregion
74:
75:
76: #region Set Dependencies
77:
78: this._dcopy = dcopy; this._dgemv = dgemv; this._drot = drot; this._xerbla = xerbla;
79:
80: #endregion
81:
82: }
83: /// <summary>
84: /// Purpose
85: /// =======
86: ///
87: /// DLAEDA computes the Z vector corresponding to the merge step in the
88: /// CURLVLth step of the merge process with TLVLS steps for the CURPBMth
89: /// problem.
90: ///
91: ///</summary>
92: /// <param name="N">
93: /// (input) INTEGER
94: /// The dimension of the symmetric tridiagonal matrix. N .GE. 0.
95: ///</param>
96: /// <param name="TLVLS">
97: /// (input) INTEGER
98: /// The total number of merging levels in the overall divide and
99: /// conquer tree.
100: ///</param>
101: /// <param name="CURLVL">
102: /// (input) INTEGER
103: /// The current level in the overall merge routine,
104: /// 0 .LE. curlvl .LE. tlvls.
105: ///</param>
106: /// <param name="CURPBM">
107: /// (input) INTEGER
108: /// The current problem in the current level in the overall
109: /// merge routine (counting from upper left to lower right).
110: ///</param>
111: /// <param name="PRMPTR">
112: /// (input) INTEGER array, dimension (N lg N)
113: /// Contains a list of pointers which indicate where in PERM a
114: /// level's permutation is stored. PRMPTR(i+1) - PRMPTR(i)
115: /// indicates the size of the permutation and incidentally the
116: /// size of the full, non-deflated problem.
117: ///</param>
118: /// <param name="PERM">
119: /// (input) INTEGER array, dimension (N lg N)
120: /// Contains the permutations (from deflation and sorting) to be
121: /// applied to each eigenblock.
122: ///</param>
123: /// <param name="GIVPTR">
124: /// (input) INTEGER array, dimension (N lg N)
125: /// Contains a list of pointers which indicate where in GIVCOL a
126: /// level's Givens rotations are stored. GIVPTR(i+1) - GIVPTR(i)
127: /// indicates the number of Givens rotations.
128: ///</param>
129: /// <param name="GIVCOL">
130: /// (input) INTEGER array, dimension (2, N lg N)
131: /// Each pair of numbers indicates a pair of columns to take place
132: /// in a Givens rotation.
133: ///</param>
134: /// <param name="GIVNUM">
135: /// (input) DOUBLE PRECISION array, dimension (2, N lg N)
136: /// Each number indicates the S value to be used in the
137: /// corresponding Givens rotation.
138: ///</param>
139: /// <param name="Q">
140: /// (input) DOUBLE PRECISION array, dimension (N**2)
141: /// Contains the square eigenblocks from previous levels, the
142: /// starting positions for blocks are given by QPTR.
143: ///</param>
144: /// <param name="QPTR">
145: /// (input) INTEGER array, dimension (N+2)
146: /// Contains a list of pointers which indicate where in Q an
147: /// eigenblock is stored. SQRT( QPTR(i+1) - QPTR(i) ) indicates
148: /// the size of the block.
149: ///</param>
150: /// <param name="Z">
151: /// (output) DOUBLE PRECISION array, dimension (N)
152: /// On output this vector contains the updating vector (the last
153: /// row of the first sub-eigenvector matrix and the first row of
154: /// the second sub-eigenvector matrix).
155: ///</param>
156: /// <param name="ZTEMP">
157: /// (workspace) DOUBLE PRECISION array, dimension (N)
158: ///</param>
159: /// <param name="INFO">
160: /// (output) INTEGER
161: /// = 0: successful exit.
162: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
163: ///</param>
164: public void Run(int N, int TLVLS, int CURLVL, int CURPBM, int[] PRMPTR, int offset_prmptr, int[] PERM, int offset_perm
165: , int[] GIVPTR, int offset_givptr, int[] GIVCOL, int offset_givcol, double[] GIVNUM, int offset_givnum, double[] Q, int offset_q, int[] QPTR, int offset_qptr, ref double[] Z, int offset_z
166: , ref double[] ZTEMP, int offset_ztemp, ref int INFO)
167: {
168:
169: #region Array Index Correction
170:
171: int o_prmptr = -1 + offset_prmptr; int o_perm = -1 + offset_perm; int o_givptr = -1 + offset_givptr;
172: int o_givcol = -3 + offset_givcol; int o_givnum = -3 + offset_givnum; int o_q = -1 + offset_q;
173: int o_qptr = -1 + offset_qptr; int o_z = -1 + offset_z; int o_ztemp = -1 + offset_ztemp;
174:
175: #endregion
176:
177:
178: #region Prolog
179:
180: // *
181: // * -- LAPACK routine (version 3.1) --
182: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
183: // * November 2006
184: // *
185: // * .. Scalar Arguments ..
186: // * ..
187: // * .. Array Arguments ..
188: // * ..
189: // *
190: // * Purpose
191: // * =======
192: // *
193: // * DLAEDA computes the Z vector corresponding to the merge step in the
194: // * CURLVLth step of the merge process with TLVLS steps for the CURPBMth
195: // * problem.
196: // *
197: // * Arguments
198: // * =========
199: // *
200: // * N (input) INTEGER
201: // * The dimension of the symmetric tridiagonal matrix. N >= 0.
202: // *
203: // * TLVLS (input) INTEGER
204: // * The total number of merging levels in the overall divide and
205: // * conquer tree.
206: // *
207: // * CURLVL (input) INTEGER
208: // * The current level in the overall merge routine,
209: // * 0 <= curlvl <= tlvls.
210: // *
211: // * CURPBM (input) INTEGER
212: // * The current problem in the current level in the overall
213: // * merge routine (counting from upper left to lower right).
214: // *
215: // * PRMPTR (input) INTEGER array, dimension (N lg N)
216: // * Contains a list of pointers which indicate where in PERM a
217: // * level's permutation is stored. PRMPTR(i+1) - PRMPTR(i)
218: // * indicates the size of the permutation and incidentally the
219: // * size of the full, non-deflated problem.
220: // *
221: // * PERM (input) INTEGER array, dimension (N lg N)
222: // * Contains the permutations (from deflation and sorting) to be
223: // * applied to each eigenblock.
224: // *
225: // * GIVPTR (input) INTEGER array, dimension (N lg N)
226: // * Contains a list of pointers which indicate where in GIVCOL a
227: // * level's Givens rotations are stored. GIVPTR(i+1) - GIVPTR(i)
228: // * indicates the number of Givens rotations.
229: // *
230: // * GIVCOL (input) INTEGER array, dimension (2, N lg N)
231: // * Each pair of numbers indicates a pair of columns to take place
232: // * in a Givens rotation.
233: // *
234: // * GIVNUM (input) DOUBLE PRECISION array, dimension (2, N lg N)
235: // * Each number indicates the S value to be used in the
236: // * corresponding Givens rotation.
237: // *
238: // * Q (input) DOUBLE PRECISION array, dimension (N**2)
239: // * Contains the square eigenblocks from previous levels, the
240: // * starting positions for blocks are given by QPTR.
241: // *
242: // * QPTR (input) INTEGER array, dimension (N+2)
243: // * Contains a list of pointers which indicate where in Q an
244: // * eigenblock is stored. SQRT( QPTR(i+1) - QPTR(i) ) indicates
245: // * the size of the block.
246: // *
247: // * Z (output) DOUBLE PRECISION array, dimension (N)
248: // * On output this vector contains the updating vector (the last
249: // * row of the first sub-eigenvector matrix and the first row of
250: // * the second sub-eigenvector matrix).
251: // *
252: // * ZTEMP (workspace) DOUBLE PRECISION array, dimension (N)
253: // *
254: // * INFO (output) INTEGER
255: // * = 0: successful exit.
256: // * < 0: if INFO = -i, the i-th argument had an illegal value.
257: // *
258: // * Further Details
259: // * ===============
260: // *
261: // * Based on contributions by
262: // * Jeff Rutter, Computer Science Division, University of California
263: // * at Berkeley, USA
264: // *
265: // * =====================================================================
266: // *
267: // * .. Parameters ..
268: // * ..
269: // * .. Local Scalars ..
270: // * ..
271: // * .. External Subroutines ..
272: // * ..
273: // * .. Intrinsic Functions ..
274: // INTRINSIC DBLE, INT, SQRT;
275: // * ..
276: // * .. Executable Statements ..
277: // *
278: // * Test the input parameters.
279: // *
280:
281: #endregion
282:
283:
284: #region Body
285:
286: INFO = 0;
287: // *
288: if (N < 0)
289: {
290: INFO = - 1;
291: }
292: if (INFO != 0)
293: {
294: this._xerbla.Run("DLAEDA", - INFO);
295: return;
296: }
297: // *
298: // * Quick return if possible
299: // *
300: if (N == 0) return;
301: // *
302: // * Determine location of first number in second half.
303: // *
304: MID = N / 2 + 1;
305: // *
306: // * Gather last/first rows of appropriate eigenblocks into center of Z
307: // *
308: PTR = 1;
309: // *
310: // * Determine location of lowest level subproblem in the full storage
311: // * scheme
312: // *
313: CURR = Convert.ToInt32(PTR + CURPBM * Math.Pow(2, CURLVL) + Math.Pow(2, CURLVL - 1) - 1);
314: // *
315: // * Determine size of these matrices. We add HALF to the value of
316: // * the SQRT in case the machine underestimates one of these square
317: // * roots.
318: // *
319: BSIZ1 = Convert.ToInt32(Math.Truncate(HALF + Math.Sqrt(Convert.ToDouble(QPTR[CURR + 1 + o_qptr] - QPTR[CURR + o_qptr]))));
320: BSIZ2 = Convert.ToInt32(Math.Truncate(HALF + Math.Sqrt(Convert.ToDouble(QPTR[CURR + 2 + o_qptr] - QPTR[CURR + 1 + o_qptr]))));
321: for (K = 1; K <= MID - BSIZ1 - 1; K++)
322: {
323: Z[K + o_z] = ZERO;
324: }
325: this._dcopy.Run(BSIZ1, Q, QPTR[CURR + o_qptr] + BSIZ1 - 1 + o_q, BSIZ1, ref Z, MID - BSIZ1 + o_z, 1);
326: this._dcopy.Run(BSIZ2, Q, QPTR[CURR + 1 + o_qptr] + o_q, BSIZ2, ref Z, MID + o_z, 1);
327: for (K = MID + BSIZ2; K <= N; K++)
328: {
329: Z[K + o_z] = ZERO;
330: }
331: // *
332: // * Loop thru remaining levels 1 -> CURLVL applying the Givens
333: // * rotations and permutation and then multiplying the center matrices
334: // * against the current Z.
335: // *
336: PTR = (int)Math.Pow(2, TLVLS) + 1;
337: for (K = 1; K <= CURLVL - 1; K++)
338: {
339: CURR = Convert.ToInt32( PTR + CURPBM * Math.Pow(2,CURLVL - K) + Math.Pow(2,CURLVL - K - 1) - 1);
340: PSIZ1 = PRMPTR[CURR + 1 + o_prmptr] - PRMPTR[CURR + o_prmptr];
341: PSIZ2 = PRMPTR[CURR + 2 + o_prmptr] - PRMPTR[CURR + 1 + o_prmptr];
342: ZPTR1 = MID - PSIZ1;
343: // *
344: // * Apply Givens at CURR and CURR+1
345: // *
346: for (I = GIVPTR[CURR + o_givptr]; I <= GIVPTR[CURR + 1 + o_givptr] - 1; I++)
347: {
348: this._drot.Run(1, ref Z, ZPTR1 + GIVCOL[1+I * 2 + o_givcol] - 1 + o_z, 1, ref Z, ZPTR1 + GIVCOL[2+I * 2 + o_givcol] - 1 + o_z, 1, GIVNUM[1+I * 2 + o_givnum]
349: , GIVNUM[2+I * 2 + o_givnum]);
350: }
351: for (I = GIVPTR[CURR + 1 + o_givptr]; I <= GIVPTR[CURR + 2 + o_givptr] - 1; I++)
352: {
353: this._drot.Run(1, ref Z, MID - 1 + GIVCOL[1+I * 2 + o_givcol] + o_z, 1, ref Z, MID - 1 + GIVCOL[2+I * 2 + o_givcol] + o_z, 1, GIVNUM[1+I * 2 + o_givnum]
354: , GIVNUM[2+I * 2 + o_givnum]);
355: }
356: PSIZ1 = PRMPTR[CURR + 1 + o_prmptr] - PRMPTR[CURR + o_prmptr];
357: PSIZ2 = PRMPTR[CURR + 2 + o_prmptr] - PRMPTR[CURR + 1 + o_prmptr];
358: for (I = 0; I <= PSIZ1 - 1; I++)
359: {
360: ZTEMP[I + 1 + o_ztemp] = Z[ZPTR1 + PERM[PRMPTR[CURR + o_prmptr] + I + o_perm] - 1 + o_z];
361: }
362: for (I = 0; I <= PSIZ2 - 1; I++)
363: {
364: ZTEMP[PSIZ1 + I + 1 + o_ztemp] = Z[MID + PERM[PRMPTR[CURR + 1 + o_prmptr] + I + o_perm] - 1 + o_z];
365: }
366: // *
367: // * Multiply Blocks at CURR and CURR+1
368: // *
369: // * Determine size of these matrices. We add HALF to the value of
370: // * the SQRT in case the machine underestimates one of these
371: // * square roots.
372: // *
373: BSIZ1 = Convert.ToInt32(Math.Truncate(HALF + Math.Sqrt(Convert.ToDouble(QPTR[CURR + 1 + o_qptr] - QPTR[CURR + o_qptr]))));
374: BSIZ2 = Convert.ToInt32(Math.Truncate(HALF + Math.Sqrt(Convert.ToDouble(QPTR[CURR + 2 + o_qptr] - QPTR[CURR + 1 + o_qptr]))));
375: if (BSIZ1 > 0)
376: {
377: this._dgemv.Run("T", BSIZ1, BSIZ1, ONE, Q, QPTR[CURR + o_qptr] + o_q, BSIZ1
378: , ZTEMP, 1 + o_ztemp, 1, ZERO, ref Z, ZPTR1 + o_z, 1);
379: }
380: this._dcopy.Run(PSIZ1 - BSIZ1, ZTEMP, BSIZ1 + 1 + o_ztemp, 1, ref Z, ZPTR1 + BSIZ1 + o_z, 1);
381: if (BSIZ2 > 0)
382: {
383: this._dgemv.Run("T", BSIZ2, BSIZ2, ONE, Q, QPTR[CURR + 1 + o_qptr] + o_q, BSIZ2
384: , ZTEMP, PSIZ1 + 1 + o_ztemp, 1, ZERO, ref Z, MID + o_z, 1);
385: }
386: this._dcopy.Run(PSIZ2 - BSIZ2, ZTEMP, PSIZ1 + BSIZ2 + 1 + o_ztemp, 1, ref Z, MID + BSIZ2 + o_z, 1);
387: // *
388: PTR = PTR + (int)Math.Pow(2, TLVLS - K);
389: }
390: // *
391: return;
392: // *
393: // * End of DLAEDA
394: // *
395:
396: #endregion
397:
398: }
399: }
400: }