~ubuntu-branches/ubuntu/lucid/igraph/lucid

« back to all changes in this revision

Viewing changes to src/arpack/dsgets.c

  • Committer: Bazaar Package Importer
  • Author(s): Mathieu Malaterre
  • Date: 2009-11-16 18:12:42 UTC
  • Revision ID: james.westby@ubuntu.com-20091116181242-mzv9p5fz9uj57xd1
Tags: upstream-0.5.3
ImportĀ upstreamĀ versionĀ 0.5.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
/*  -- translated by f2c (version 20050501).
 
2
   You must link the resulting object file with libf2c:
 
3
        on Microsoft Windows system, link with libf2c.lib;
 
4
        on Linux or Unix systems, link with .../path/to/libf2c.a -lm
 
5
        or, if you install libf2c.a in a standard place, with -lf2c -lm
 
6
        -- in that order, at the end of the command line, as in
 
7
                cc *.o -lf2c -lm
 
8
        Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
 
9
 
 
10
                http://www.netlib.org/f2c/libf2c.zip
 
11
*/
 
12
 
 
13
#include "config.h"
 
14
#include "arpack_internal.h"
 
15
#include "f2c.h"
 
16
 
 
17
/* Common Block Declarations */
 
18
 
 
19
static struct {
 
20
    integer logfil, ndigit, mgetv0, msaupd, msaup2, msaitr, mseigt, msapps, 
 
21
            msgets, mseupd, mnaupd, mnaup2, mnaitr, mneigh, mnapps, mngets, 
 
22
            mneupd, mcaupd, mcaup2, mcaitr, mceigh, mcapps, mcgets, mceupd;
 
23
} debug_;
 
24
 
 
25
#define debug_1 debug_
 
26
 
 
27
static struct {
 
28
    integer nopx, nbx, nrorth, nitref, nrstrt;
 
29
    real tsaupd, tsaup2, tsaitr, tseigt, tsgets, tsapps, tsconv, tnaupd, 
 
30
            tnaup2, tnaitr, tneigh, tngets, tnapps, tnconv, tcaupd, tcaup2, 
 
31
            tcaitr, tceigh, tcgets, tcapps, tcconv, tmvopx, tmvbx, tgetv0, 
 
32
            titref, trvec;
 
33
} timing_;
 
34
 
 
35
#define timing_1 timing_
 
36
 
 
37
/* Table of constant values */
 
38
 
 
39
static logical c_true = TRUE_;
 
40
static integer c__1 = 1;
 
41
 
 
42
/* ----------------------------------------------------------------------- */
 
43
/* \BeginDoc */
 
44
 
 
45
/* \Name: dsgets */
 
46
 
 
47
/* \Description: */
 
48
/*  Given the eigenvalues of the symmetric tridiagonal matrix H, */
 
49
/*  computes the NP shifts AMU that are zeros of the polynomial of */
 
50
/*  degree NP which filters out components of the unwanted eigenvectors */
 
51
/*  corresponding to the AMU's based on some given criteria. */
 
52
 
 
53
/*  NOTE: This is called even in the case of user specified shifts in */
 
54
/*  order to sort the eigenvalues, and error bounds of H for later use. */
 
55
 
 
56
/* \Usage: */
 
57
/*  call dsgets */
 
58
/*     ( ISHIFT, WHICH, KEV, NP, RITZ, BOUNDS, SHIFTS ) */
 
59
 
 
60
/* \Arguments */
 
61
/*  ISHIFT  Integer.  (INPUT) */
 
62
/*          Method for selecting the implicit shifts at each iteration. */
 
63
/*          ISHIFT = 0: user specified shifts */
 
64
/*          ISHIFT = 1: exact shift with respect to the matrix H. */
 
65
 
 
66
/*  WHICH   Character*2.  (INPUT) */
 
67
/*          Shift selection criteria. */
 
68
/*          'LM' -> KEV eigenvalues of largest magnitude are retained. */
 
69
/*          'SM' -> KEV eigenvalues of smallest magnitude are retained. */
 
70
/*          'LA' -> KEV eigenvalues of largest value are retained. */
 
71
/*          'SA' -> KEV eigenvalues of smallest value are retained. */
 
72
/*          'BE' -> KEV eigenvalues, half from each end of the spectrum. */
 
73
/*                  If KEV is odd, compute one more from the high end. */
 
74
 
 
75
/*  KEV      Integer.  (INPUT) */
 
76
/*          KEV+NP is the size of the matrix H. */
 
77
 
 
78
/*  NP      Integer.  (INPUT) */
 
79
/*          Number of implicit shifts to be computed. */
 
80
 
 
81
/*  RITZ    Double precision array of length KEV+NP.  (INPUT/OUTPUT) */
 
82
/*          On INPUT, RITZ contains the eigenvalues of H. */
 
83
/*          On OUTPUT, RITZ are sorted so that the unwanted eigenvalues */
 
84
/*          are in the first NP locations and the wanted part is in */
 
85
/*          the last KEV locations.  When exact shifts are selected, the */
 
86
/*          unwanted part corresponds to the shifts to be applied. */
 
87
 
 
88
/*  BOUNDS  Double precision array of length KEV+NP.  (INPUT/OUTPUT) */
 
89
/*          Error bounds corresponding to the ordering in RITZ. */
 
90
 
 
91
/*  SHIFTS  Double precision array of length NP.  (INPUT/OUTPUT) */
 
92
/*          On INPUT:  contains the user specified shifts if ISHIFT = 0. */
 
93
/*          On OUTPUT: contains the shifts sorted into decreasing order */
 
94
/*          of magnitude with respect to the Ritz estimates contained in */
 
95
/*          BOUNDS. If ISHIFT = 0, SHIFTS is not modified on exit. */
 
96
 
 
97
/* \EndDoc */
 
98
 
 
99
/* ----------------------------------------------------------------------- */
 
100
 
 
101
/* \BeginLib */
 
102
 
 
103
/* \Local variables: */
 
104
/*     xxxxxx  real */
 
105
 
 
106
/* \Routines called: */
 
107
/*     dsortr  ARPACK utility sorting routine. */
 
108
/*     ivout   ARPACK utility routine that prints integers. */
 
109
/*     second  ARPACK utility routine for timing. */
 
110
/*     dvout   ARPACK utility routine that prints vectors. */
 
111
/*     dcopy   Level 1 BLAS that copies one vector to another. */
 
112
/*     dswap   Level 1 BLAS that swaps the contents of two vectors. */
 
113
 
 
114
/* \Author */
 
115
/*     Danny Sorensen               Phuong Vu */
 
116
/*     Richard Lehoucq              CRPC / Rice University */
 
117
/*     Dept. of Computational &     Houston, Texas */
 
118
/*     Applied Mathematics */
 
119
/*     Rice University */
 
120
/*     Houston, Texas */
 
121
 
 
122
/* \Revision history: */
 
123
/*     xx/xx/93: Version ' 2.1' */
 
124
 
 
125
/* \SCCS Information: @(#) */
 
126
/* FILE: sgets.F   SID: 2.4   DATE OF SID: 4/19/96   RELEASE: 2 */
 
127
 
 
128
/* \Remarks */
 
129
 
 
130
/* \EndLib */
 
131
 
 
132
/* ----------------------------------------------------------------------- */
 
133
 
 
134
/* Subroutine */ int igraphdsgets_(integer *ishift, char *which, integer *kev,
 
135
         integer *np, doublereal *ritz, doublereal *bounds, doublereal *
 
136
        shifts)
 
137
{
 
138
    /* System generated locals */
 
139
    integer i__1;
 
140
 
 
141
    /* Builtin functions */
 
142
    integer igraphs_cmp(char *, char *, ftnlen, ftnlen);
 
143
 
 
144
    /* Local variables */
 
145
    static real t0, t1;
 
146
    extern /* Subroutine */ int igraphdswap_(integer *, doublereal *, integer 
 
147
            *, doublereal *, integer *), igraphdcopy_(integer *, doublereal *,
 
148
             integer *, doublereal *, integer *), igraphdvout_(integer *, 
 
149
            integer *, doublereal *, integer *, char *), igraphivout_(
 
150
            integer *, integer *, integer *, integer *, char *), 
 
151
            igraphsecond_(real *);
 
152
    static integer kevd2;
 
153
    extern /* Subroutine */ int igraphdsortr_(char *, logical *, integer *, 
 
154
            doublereal *, doublereal *);
 
155
    static integer msglvl;
 
156
 
 
157
 
 
158
/*     %----------------------------------------------------% */
 
159
/*     | Include files for debugging and timing information | */
 
160
/*     %----------------------------------------------------% */
 
161
 
 
162
 
 
163
/* \SCCS Information: @(#) */
 
164
/* FILE: debug.h   SID: 2.3   DATE OF SID: 11/16/95   RELEASE: 2 */
 
165
 
 
166
/*     %---------------------------------% */
 
167
/*     | See debug.doc for documentation | */
 
168
/*     %---------------------------------% */
 
169
 
 
170
/*     %------------------% */
 
171
/*     | Scalar Arguments | */
 
172
/*     %------------------% */
 
173
 
 
174
/*     %--------------------------------% */
 
175
/*     | See stat.doc for documentation | */
 
176
/*     %--------------------------------% */
 
177
 
 
178
/* \SCCS Information: @(#) */
 
179
/* FILE: stat.h   SID: 2.2   DATE OF SID: 11/16/95   RELEASE: 2 */
 
180
 
 
181
 
 
182
 
 
183
/*     %-----------------% */
 
184
/*     | Array Arguments | */
 
185
/*     %-----------------% */
 
186
 
 
187
 
 
188
/*     %------------% */
 
189
/*     | Parameters | */
 
190
/*     %------------% */
 
191
 
 
192
 
 
193
/*     %---------------% */
 
194
/*     | Local Scalars | */
 
195
/*     %---------------% */
 
196
 
 
197
 
 
198
/*     %----------------------% */
 
199
/*     | External Subroutines | */
 
200
/*     %----------------------% */
 
201
 
 
202
 
 
203
/*     %---------------------% */
 
204
/*     | Intrinsic Functions | */
 
205
/*     %---------------------% */
 
206
 
 
207
 
 
208
/*     %-----------------------% */
 
209
/*     | Executable Statements | */
 
210
/*     %-----------------------% */
 
211
 
 
212
/*     %-------------------------------% */
 
213
/*     | Initialize timing statistics  | */
 
214
/*     | & message level for debugging | */
 
215
/*     %-------------------------------% */
 
216
 
 
217
    /* Parameter adjustments */
 
218
    --shifts;
 
219
    --bounds;
 
220
    --ritz;
 
221
 
 
222
    /* Function Body */
 
223
    igraphsecond_(&t0);
 
224
    msglvl = debug_1.msgets;
 
225
 
 
226
    if (igraphs_cmp(which, "BE", (ftnlen)2, (ftnlen)2) == 0) {
 
227
 
 
228
/*        %-----------------------------------------------------% */
 
229
/*        | Both ends of the spectrum are requested.            | */
 
230
/*        | Sort the eigenvalues into algebraically increasing  | */
 
231
/*        | order first then swap high end of the spectrum next | */
 
232
/*        | to low end in appropriate locations.                | */
 
233
/*        | NOTE: when np < floor(kev/2) be careful not to swap | */
 
234
/*        | overlapping locations.                              | */
 
235
/*        %-----------------------------------------------------% */
 
236
 
 
237
        i__1 = *kev + *np;
 
238
        igraphdsortr_("LA", &c_true, &i__1, &ritz[1], &bounds[1]);
 
239
        kevd2 = *kev / 2;
 
240
        if (*kev > 1) {
 
241
            i__1 = min(kevd2,*np);
 
242
            igraphdswap_(&i__1, &ritz[1], &c__1, &ritz[max(kevd2,*np) + 1], &
 
243
                    c__1);
 
244
            i__1 = min(kevd2,*np);
 
245
            igraphdswap_(&i__1, &bounds[1], &c__1, &bounds[max(kevd2,*np) + 1]
 
246
                    , &c__1);
 
247
        }
 
248
 
 
249
    } else {
 
250
 
 
251
/*        %----------------------------------------------------% */
 
252
/*        | LM, SM, LA, SA case.                               | */
 
253
/*        | Sort the eigenvalues of H into the desired order   | */
 
254
/*        | and apply the resulting order to BOUNDS.           | */
 
255
/*        | The eigenvalues are sorted so that the wanted part | */
 
256
/*        | are always in the last KEV locations.               | */
 
257
/*        %----------------------------------------------------% */
 
258
 
 
259
        i__1 = *kev + *np;
 
260
        igraphdsortr_(which, &c_true, &i__1, &ritz[1], &bounds[1]);
 
261
    }
 
262
 
 
263
    if (*ishift == 1 && *np > 0) {
 
264
 
 
265
/*        %-------------------------------------------------------% */
 
266
/*        | Sort the unwanted Ritz values used as shifts so that  | */
 
267
/*        | the ones with largest Ritz estimates are first.       | */
 
268
/*        | This will tend to minimize the effects of the         | */
 
269
/*        | forward instability of the iteration when the shifts  | */
 
270
/*        | are applied in subroutine dsapps.                     | */
 
271
/*        %-------------------------------------------------------% */
 
272
 
 
273
        igraphdsortr_("SM", &c_true, np, &bounds[1], &ritz[1]);
 
274
        igraphdcopy_(np, &ritz[1], &c__1, &shifts[1], &c__1);
 
275
    }
 
276
 
 
277
    igraphsecond_(&t1);
 
278
    timing_1.tsgets += t1 - t0;
 
279
 
 
280
    if (msglvl > 0) {
 
281
        igraphivout_(&debug_1.logfil, &c__1, kev, &debug_1.ndigit, "_sgets: "
 
282
                "KEV is");
 
283
        igraphivout_(&debug_1.logfil, &c__1, np, &debug_1.ndigit, "_sgets: N"
 
284
                "P is");
 
285
        i__1 = *kev + *np;
 
286
        igraphdvout_(&debug_1.logfil, &i__1, &ritz[1], &debug_1.ndigit, "_sg"
 
287
                "ets: Eigenvalues of current H matrix");
 
288
        i__1 = *kev + *np;
 
289
        igraphdvout_(&debug_1.logfil, &i__1, &bounds[1], &debug_1.ndigit, 
 
290
                "_sgets: Associated Ritz estimates");
 
291
    }
 
292
 
 
293
    return 0;
 
294
 
 
295
/*     %---------------% */
 
296
/*     | End of dsgets | */
 
297
/*     %---------------% */
 
298
 
 
299
} /* igraphdsgets_ */
 
300