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.0) --
21: /// Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
22: /// Courant Institute, Argonne National Lab, and Rice University
23: /// September 30, 1994
24: /// Purpose
25: /// =======
26: ///
27: /// DDISNA computes the reciprocal condition numbers for the eigenvectors
28: /// of a real symmetric or complex Hermitian matrix or for the left or
29: /// right singular vectors of a general m-by-n matrix. The reciprocal
30: /// condition number is the 'gap' between the corresponding eigenvalue or
31: /// singular value and the nearest other one.
32: ///
33: /// The bound on the error, measured by angle in radians, in the I-th
34: /// computed vector is given by
35: ///
36: /// DLAMCH( 'E' ) * ( ANORM / SEP( I ) )
37: ///
38: /// where ANORM = 2-norm(A) = max( abs( D(j) ) ). SEP(I) is not allowed
39: /// to be smaller than DLAMCH( 'E' )*ANORM in order to limit the size of
40: /// the error bound.
41: ///
42: /// DDISNA may also be used to compute error bounds for eigenvectors of
43: /// the generalized symmetric definite eigenproblem.
44: ///
45: ///</summary>
46: public class DDISNA
47: {
48:
49:
50: #region Dependencies
51:
52: LSAME _lsame; DLAMCH _dlamch; XERBLA _xerbla;
53:
54: #endregion
55:
56:
57: #region Fields
58:
59: const double ZERO = 0.0E+0; bool DECR = false; bool EIGEN = false; bool INCR = false; bool LEFT = false;
60: bool RIGHT = false;bool SING = false; int I = 0; int K = 0; double ANORM = 0; double EPS = 0; double NEWGAP = 0;
61: double OLDGAP = 0;double SAFMIN = 0; double THRESH = 0;
62:
63: #endregion
64:
65: public DDISNA(LSAME lsame, DLAMCH dlamch, XERBLA xerbla)
66: {
67:
68:
69: #region Set Dependencies
70:
71: this._lsame = lsame; this._dlamch = dlamch; this._xerbla = xerbla;
72:
73: #endregion
74:
75: }
76:
77: public DDISNA()
78: {
79:
80:
81: #region Dependencies (Initialization)
82:
83: LSAME lsame = new LSAME();
84: DLAMC3 dlamc3 = new DLAMC3();
85: XERBLA xerbla = new XERBLA();
86: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
87: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
88: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
89: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
90: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
91:
92: #endregion
93:
94:
95: #region Set Dependencies
96:
97: this._lsame = lsame; this._dlamch = dlamch; this._xerbla = xerbla;
98:
99: #endregion
100:
101: }
102: /// <summary>
103: /// Purpose
104: /// =======
105: ///
106: /// DDISNA computes the reciprocal condition numbers for the eigenvectors
107: /// of a real symmetric or complex Hermitian matrix or for the left or
108: /// right singular vectors of a general m-by-n matrix. The reciprocal
109: /// condition number is the 'gap' between the corresponding eigenvalue or
110: /// singular value and the nearest other one.
111: ///
112: /// The bound on the error, measured by angle in radians, in the I-th
113: /// computed vector is given by
114: ///
115: /// DLAMCH( 'E' ) * ( ANORM / SEP( I ) )
116: ///
117: /// where ANORM = 2-norm(A) = max( abs( D(j) ) ). SEP(I) is not allowed
118: /// to be smaller than DLAMCH( 'E' )*ANORM in order to limit the size of
119: /// the error bound.
120: ///
121: /// DDISNA may also be used to compute error bounds for eigenvectors of
122: /// the generalized symmetric definite eigenproblem.
123: ///
124: ///</summary>
125: /// <param name="JOB">
126: /// (input) CHARACTER*1
127: /// Specifies for which problem the reciprocal condition numbers
128: /// should be computed:
129: /// = 'E': the eigenvectors of a symmetric/Hermitian matrix;
130: /// = 'L': the left singular vectors of a general matrix;
131: /// = 'R': the right singular vectors of a general matrix.
132: ///</param>
133: /// <param name="M">
134: /// (input) INTEGER
135: /// The number of rows of the matrix. M .GE. 0.
136: ///</param>
137: /// <param name="N">
138: /// (input) INTEGER
139: /// If JOB = 'L' or 'R', the number of columns of the matrix,
140: /// in which case N .GE. 0. Ignored if JOB = 'E'.
141: ///</param>
142: /// <param name="D">
143: /// (input) DOUBLE PRECISION array, dimension (M) if JOB = 'E'
144: /// dimension (min(M,N)) if JOB = 'L' or 'R'
145: /// The eigenvalues (if JOB = 'E') or singular values (if JOB =
146: /// 'L' or 'R') of the matrix, in either increasing or decreasing
147: /// order. If singular values, they must be non-negative.
148: ///</param>
149: /// <param name="SEP">
150: /// (output) DOUBLE PRECISION array, dimension (M) if JOB = 'E'
151: /// dimension (min(M,N)) if JOB = 'L' or 'R'
152: /// The reciprocal condition numbers of the vectors.
153: ///</param>
154: /// <param name="INFO">
155: /// (output) INTEGER
156: /// = 0: successful exit.
157: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
158: ///</param>
159: public void Run(string JOB, int M, int N, double[] D, int offset_d, ref double[] SEP, int offset_sep, ref int INFO)
160: {
161:
162: #region Array Index Correction
163:
164: int o_d = -1 + offset_d; int o_sep = -1 + offset_sep;
165:
166: #endregion
167:
168:
169: #region Strings
170:
171: JOB = JOB.Substring(0, 1);
172:
173: #endregion
174:
175:
176: #region Prolog
177:
178: // *
179: // * -- LAPACK routine (version 3.0) --
180: // * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
181: // * Courant Institute, Argonne National Lab, and Rice University
182: // * September 30, 1994
183: // *
184: // * .. Scalar Arguments ..
185: // * ..
186: // * .. Array Arguments ..
187: // * ..
188: // *
189: // * Purpose
190: // * =======
191: // *
192: // * DDISNA computes the reciprocal condition numbers for the eigenvectors
193: // * of a real symmetric or complex Hermitian matrix or for the left or
194: // * right singular vectors of a general m-by-n matrix. The reciprocal
195: // * condition number is the 'gap' between the corresponding eigenvalue or
196: // * singular value and the nearest other one.
197: // *
198: // * The bound on the error, measured by angle in radians, in the I-th
199: // * computed vector is given by
200: // *
201: // * DLAMCH( 'E' ) * ( ANORM / SEP( I ) )
202: // *
203: // * where ANORM = 2-norm(A) = max( abs( D(j) ) ). SEP(I) is not allowed
204: // * to be smaller than DLAMCH( 'E' )*ANORM in order to limit the size of
205: // * the error bound.
206: // *
207: // * DDISNA may also be used to compute error bounds for eigenvectors of
208: // * the generalized symmetric definite eigenproblem.
209: // *
210: // * Arguments
211: // * =========
212: // *
213: // * JOB (input) CHARACTER*1
214: // * Specifies for which problem the reciprocal condition numbers
215: // * should be computed:
216: // * = 'E': the eigenvectors of a symmetric/Hermitian matrix;
217: // * = 'L': the left singular vectors of a general matrix;
218: // * = 'R': the right singular vectors of a general matrix.
219: // *
220: // * M (input) INTEGER
221: // * The number of rows of the matrix. M >= 0.
222: // *
223: // * N (input) INTEGER
224: // * If JOB = 'L' or 'R', the number of columns of the matrix,
225: // * in which case N >= 0. Ignored if JOB = 'E'.
226: // *
227: // * D (input) DOUBLE PRECISION array, dimension (M) if JOB = 'E'
228: // * dimension (min(M,N)) if JOB = 'L' or 'R'
229: // * The eigenvalues (if JOB = 'E') or singular values (if JOB =
230: // * 'L' or 'R') of the matrix, in either increasing or decreasing
231: // * order. If singular values, they must be non-negative.
232: // *
233: // * SEP (output) DOUBLE PRECISION array, dimension (M) if JOB = 'E'
234: // * dimension (min(M,N)) if JOB = 'L' or 'R'
235: // * The reciprocal condition numbers of the vectors.
236: // *
237: // * INFO (output) INTEGER
238: // * = 0: successful exit.
239: // * < 0: if INFO = -i, the i-th argument had an illegal value.
240: // *
241: // * =====================================================================
242: // *
243: // * .. Parameters ..
244: // * ..
245: // * .. Local Scalars ..
246: // * ..
247: // * .. External Functions ..
248: // * ..
249: // * .. Intrinsic Functions ..
250: // INTRINSIC ABS, MAX, MIN;
251: // * ..
252: // * .. External Subroutines ..
253: // * ..
254: // * .. Executable Statements ..
255: // *
256: // * Test the input arguments
257: // *
258:
259: #endregion
260:
261:
262: #region Body
263:
264: INFO = 0;
265: EIGEN = this._lsame.Run(JOB, "E");
266: LEFT = this._lsame.Run(JOB, "L");
267: RIGHT = this._lsame.Run(JOB, "R");
268: SING = LEFT || RIGHT;
269: if (EIGEN)
270: {
271: K = M;
272: }
273: else
274: {
275: if (SING)
276: {
277: K = Math.Min(M, N);
278: }
279: }
280: if (!EIGEN && !SING)
281: {
282: INFO = - 1;
283: }
284: else
285: {
286: if (M < 0)
287: {
288: INFO = - 2;
289: }
290: else
291: {
292: if (K < 0)
293: {
294: INFO = - 3;
295: }
296: else
297: {
298: INCR = true;
299: DECR = true;
300: for (I = 1; I <= K - 1; I++)
301: {
302: if (INCR) INCR = INCR && D[I + o_d] <= D[I + 1 + o_d];
303: if (DECR) DECR = DECR && D[I + o_d] >= D[I + 1 + o_d];
304: }
305: if (SING && K > 0)
306: {
307: if (INCR) INCR = INCR && ZERO <= D[1 + o_d];
308: if (DECR) DECR = DECR && D[K + o_d] >= ZERO;
309: }
310: if (!(INCR || DECR)) INFO = - 4;
311: }
312: }
313: }
314: if (INFO != 0)
315: {
316: this._xerbla.Run("DDISNA", - INFO);
317: return;
318: }
319: // *
320: // * Quick return if possible
321: // *
322: if (K == 0) return;
323: // *
324: // * Compute reciprocal condition numbers
325: // *
326: if (K == 1)
327: {
328: SEP[1 + o_sep] = this._dlamch.Run("O");
329: }
330: else
331: {
332: OLDGAP = Math.Abs(D[2 + o_d] - D[1 + o_d]);
333: SEP[1 + o_sep] = OLDGAP;
334: for (I = 2; I <= K - 1; I++)
335: {
336: NEWGAP = Math.Abs(D[I + 1 + o_d] - D[I + o_d]);
337: SEP[I + o_sep] = Math.Min(OLDGAP, NEWGAP);
338: OLDGAP = NEWGAP;
339: }
340: SEP[K + o_sep] = OLDGAP;
341: }
342: if (SING)
343: {
344: if ((LEFT && M > N) || (RIGHT && M < N))
345: {
346: if (INCR) SEP[1 + o_sep] = Math.Min(SEP[1 + o_sep], D[1 + o_d]);
347: if (DECR) SEP[K + o_sep] = Math.Min(SEP[K + o_sep], D[K + o_d]);
348: }
349: }
350: // *
351: // * Ensure that reciprocal condition numbers are not less than
352: // * threshold, in order to limit the size of the error bound
353: // *
354: EPS = this._dlamch.Run("E");
355: SAFMIN = this._dlamch.Run("S");
356: ANORM = Math.Max(Math.Abs(D[1 + o_d]), Math.Abs(D[K + o_d]));
357: if (ANORM == ZERO)
358: {
359: THRESH = EPS;
360: }
361: else
362: {
363: THRESH = Math.Max(EPS * ANORM, SAFMIN);
364: }
365: for (I = 1; I <= K; I++)
366: {
367: SEP[I + o_sep] = Math.Max(SEP[I + o_sep], THRESH);
368: }
369: // *
370: return;
371: // *
372: // * End of DDISNA
373: // *
374:
375: #endregion
376:
377: }
378: }
379: }