tensor-0.1.0
 All Data Structures Namespaces Functions Variables Typedefs Enumerations Enumerator Groups Pages
saupp.h
1 // -*- mode: c++; fill-column: 80; c-basic-offset: 2; indent-tabs-mode: nil -*-
2 /*
3  Copyright (c) 2010 Juan Jose Garcia Ripoll
4 
5  Tensor is free software; you can redistribute it and/or modify it
6  under the terms of the GNU Library General Public License as published
7  by the Free Software Foundation; either version 2 of the License, or
8  (at your option) any later version.
9 
10  This program is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU Library General Public License for more details.
14 
15  You should have received a copy of the GNU General Public License along
16  with this program; if not, write to the Free Software Foundation, Inc.,
17  51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
18 */
19 /*
20  ARPACK++ v1.0 8/1/1997
21  c++ interface to ARPACK code.
22 
23  MODULE saupp.h.
24  Interface to ARPACK subroutines dsaupd and ssaupd.
25 
26  ARPACK Authors
27  Richard Lehoucq
28  Danny Sorensen
29  Chao Yang
30  Dept. of Computational & Applied Mathematics
31  Rice University
32  Houston, Texas
33 */
34 
35 #ifndef SAUPP_H
36 #define SAUPP_H
37 
38 #include "arpackf.h"
39 
40 inline void saupp(integer& ido, char bmat, integer n, const char* which, integer nev,
41  double& tol, double resid[], integer ncv, double V[],
42  integer ldv, integer iparam[], integer ipntr[], double workd[],
43  double workl[], integer lworkl, integer& info)
44 
45 /*
46  c++ version of ARPACK routine dsaupd that implements a variant of
47  the Lanczos method. This method has been designed to compute
48  approximations to a few eigenpairs of a linear operator OP that is
49  real and symmetric with respect to a real positive semi-definite
50  symmetric matrix B, i.e.
51 
52  B*OP = (OP')*B.
53 
54  where A' denotes transpose of A. In the standard eigenproblem B is
55  the identity matrix. Another way to express this condition is
56 
57  < x,OPy > = < OPx,y > where < z,w > = z'Bw.
58 
59  The computed approximate eigenvalues are called Ritz values and
60  the corresponding approximate eigenvectors are called Ritz vectors.
61 
62  saupp is usually called iteratively to solve one of the
63  following problems:
64 
65  Mode 1: A*x = lambda*x, A symmetric
66  ===> OP = A and B = I.
67 
68  Mode 2: A*x = lambda*M*x, A symmetric, M symmetric positive definite
69  ===> OP = inv[M]*A and B = M.
70  ===> (If M can be factored see remark 3 below)
71 
72  Mode 3: K*x = lambda*M*x, K symmetric, M symmetric positive semi-definite
73  ===> OP = (inv[K - sigma*M])*M and B = M.
74  ===> Shift-and-Invert mode
75 
76  Mode 4: K*x = lambda*KG*x, K symmetric positive semi-definite,
77  KG symmetric indefinite
78  ===> OP = (inv[K - sigma*KG])*K and B = K.
79  ===> Buckling mode
80 
81  Mode 5: A*x = lambda*M*x, A symmetric, M symmetric positive semi-definite
82  ===> OP = inv[A - sigma*M]*[A + sigma*M] and B = M.
83  ===> Cayley transformed mode
84 
85  NOTE: The action of w <- inv[A - sigma*M]*v or w <- inv[M]*v should be
86  accomplished either by a direct method using a sparse matrix
87  factorization and solving
88 
89  [A - sigma*M]*w = v or M*w = v,
90 
91  or through an iterative method for solving these systems. If an
92  iterative method is used, the convergence test must be more
93  stringent than the accuracy requirements for the eigenvalue
94  approximations.
95 
96  Parameters:
97 
98  ido (Input / Output) Reverse communication flag. ido must be
99  zero on the first call to saupp. ido will be set
100  internally to indicate the type of operation to be
101  performed. Control is then given back to the calling
102  routine which has the responsibility to carry out the
103  requested operation and call saupp with the result. The
104  operand is given in workd[ipntr[0]], the result must be
105  put in workd[ipntr[2]]. (If Mode = 2 see remark 5 below).
106  ido = 0: first call to the reverse communication interface.
107  ido = -1: compute Y = OP * X where
108  ipntr[0] is the pointer into workd for X,
109  ipntr[2] is the pointer into workd for Y.
110  This is for the initialization phase to force the
111  starting vector into the range of OP.
112  ido = 1: compute Y = OP * X where
113  ipntr[0] is the pointer into workd for X,
114  ipntr[2] is the pointer into workd for Y.
115  In mode 3,4 and 5, the vector B * X is already
116  available in workd[ipntr[3]]. It does not
117  need to be recomputed in forming OP * X.
118  ido = 2: compute Y = B * X where
119  ipntr[0] is the pointer into workd for X,
120  ipntr[2] is the pointer into workd for Y.
121  ido = 3: compute the iparam[8] shifts where
122  ipntr[11] is the pointer into workl for
123  placing the shifts. See remark 6 below.
124  ido = 99: done.
125  bmat (Input) bmat specifies the type of the matrix B that defines
126  the semi-inner product for the operator OP.
127  bmat = 'I' -> standard eigenvalue problem A*x = lambda*x;
128  bmat = 'G' -> generalized eigenvalue problem A*x = lambda*B*x.
129  n (Input) Dimension of the eigenproblem.
130  nev (Input) Number of eigenvalues to be computed. 0 < nev < n.
131  which (Input) Specify which of the Ritz values of OP to compute.
132  'LA' - compute the nev largest (algebraic) eigenvalues.
133  'SA' - compute the nev smallest (algebraic) eigenvalues.
134  'LM' - compute the nev largest (in magnitude) eigenvalues.
135  'SM' - compute the nev smallest (in magnitude) eigenvalues.
136  'BE' - compute nev eigenvalues, half from each end of the
137  spectrum. When NEV is odd, compute one more from the
138  high end than from the low end.
139  (see remark 1 below)
140  tol (Input) Stopping criterion: the relative accuracy of the
141  Ritz value is considered acceptable if BOUNDS[i] <=
142  tol*abs(RITZ[i]). If tol<=0.0 is passed, the machine
143  precision as computed by the LAPACK auxiliary subroutine
144  _LAMCH is used.
145  resid (Input / Output) Array of length n.
146  On input:
147  If info==0, a random initial residual vector is used.
148  If info!=0, resid contains the initial residual vector,
149  possibly from a previous run.
150  On output:
151  resid contains the final residual vector.
152  ncv (Input) Number of Lanczos vectors that are generated at each
153  iteration. After the startup phase in which nev Lanczos
154  vectors are generated, the algorithm generates ncv-nev
155  Lanczos vectors at each subsequent update iteration. Most of
156  the cost in generating each Lanczos vector is in the
157  matrix-vector product OP*x. (See remark 4 below).
158  V (Output) Double precision array of length ncv*n+1. V contains
159  the ncv Lanczos basis vectors. The first element V[0] is never
160  referenced.
161  ldv (Input) Dimension of the basis vectors contianed in V. This
162  parameter MUST be set to n.
163  iparam (Input / Output) Array of length 12.
164  iparam[0] = ISHIFT: method for selecting the implicit shifts.
165  The shifts selected at each iteration are used to restart
166  the Arnoldi iteration in an implicit fashion.
167  -------------------------------------------------------------
168  ISHIFT = 0: the shifts are provided by the user via
169  reverse communication. The NCV eigenvalues of
170  the current tridiagonal matrix T are returned in
171  the part of workl array corresponding to RITZ.
172  See remark 6 below.
173  ISHIFT = 1: exact shifts with respect to the reduced
174  tridiagonal matrix T. This is equivalent to
175  restarting the iteration with a starting vector
176  that is a linear combination of Ritz vectors
177  associated with the "wanted" Ritz values.
178  -------------------------------------------------------------
179  iparam[2] is no longer referenced.
180  iparam[3] = MXITER
181  On INPUT: maximum number of Arnoldi update iterations allowed.
182  On OUTPUT: actual number of Arnoldi update iterations taken.
183  iparam[4] = NB: blocksize to be used in the recurrence.
184  The code currently works only for NB = 1.
185  iparam[5] = NCONV: number of "converged" Ritz values.
186  This represents the number of Ritz values that satisfy
187  the convergence criterion.
188  iparam[6] is no longer referenced.
189  iparam[7] = MODE. On INPUT determines what type of
190  eigenproblem is being solved. Must be 1,2,3,4,5.
191  iparam[8] = NP. When ido = 3 and the user provides shifts
192  through reverse communication (iparam[0]=0), saupp returns
193  NP, the number of shifts the user is to provide.
194  0 < NP <=ncv-nev. See Remark 6 below.
195  iparam[9] = total number of OP*x operations.
196  iparam[10] = total number of B*x operations if bmat='G'.
197  iparam[11] = total number of steps of re-orthogonalization.
198  ipntr (Output) Array of length 12. Pointer to mark the starting
199  locations in the workd and workl arrays for matrices/vectors
200  used by the Lanczos iteration.
201  ipntr[0] : pointer to the current operand vector X in workd.
202  ipntr[2] : pointer to the current result vector Y in workd.
203  ipntr[3] : pointer to the vector B * X in workd when used in
204  the shift-and-invert mode.
205  ipntr[4] : pointer to the next available location in workl
206  that is untouched by the program.
207  ipntr[5] : pointer to the ncv by 2 tridiagonal matrix T in
208  workl.
209  ipntr[6] : pointer to the ncv RITZ values array in workl.
210  ipntr[7] : pointer to the Ritz estimates in array workl
211  associated with the Ritz values located in RITZ
212  in workl.
213  ipntr[11]: pointer to the np shifts in workl. See Remark 6.
214  Note: ipntr[8:10] is only referenced by seupp. See Remark 2.
215  ipntr[8] : pointer to the ncv RITZ values of the original
216  system.
217  ipntr[9] : pointer to the ncv corresponding error bounds.
218  ipntr[10]: pointer to the ncv by ncv matrix of eigenvectors
219  of the tridiagonal matrix T. Only referenced by
220  seupp if RVEC = TRUE. See Remarks.
221  workd (Input / Output) Array of length 3*N+1.
222  Distributed array to be used in the basic Arnoldi iteration
223  for reverse communication. The user should not use workd as
224  temporary workspace during the iteration. Upon termination
225  workd[1:n] contains B*resid[1:n]. If the Ritz vectors are
226  desired subroutine seupp uses this output.
227  workl (Output) Array of length lworkl+1. Private (replicated) array
228  on each PE or array allocated on the front end.
229  lworkl (Input) lworkl must be at least ncv*(ncv+8).
230  info (Input / Output) On input, if info = 0, a randomly initial
231  residual vector is used, otherwise resid contains the initial
232  residual vector, possibly from a previous run.
233  On output, info works as a error flag:
234  = 0 : Normal exit.
235  = 1 : Maximum number of iterations taken. All possible
236  eigenvalues of OP has been found. iparam[5]
237  returns the number of wanted converged Ritz values.
238  = 3 : No shifts could be applied during a cycle of the
239  Implicitly restarted Arnoldi iteration. One
240  possibility is to increase the size of NCV relative
241  to nev. See remark 4 below.
242  = -1 : n must be positive.
243  = -2 : nev must be positive.
244  = -3 : ncv must satisfy nev < ncv <= n.
245  = -4 : The maximum number of Arnoldi update iterations allowed
246  must be greater than zero.
247  = -5 : which must be one of 'LM', 'SM', 'LA', 'SA' or 'BE'.
248  = -6 : bmat must be one of 'I' or 'G'.
249  = -7 : Length of private work array workl is not sufficient.
250  = -8 : Error return from trid. eigenvalue calculation;
251  Informational error from LAPACK routine dsteqr.
252  = -9 : Starting vector is zero.
253  = -10 : iparam[7] must be 1,2,3,4,5.
254  = -11 : iparam[7] = 1 and bmat = 'G' are incompatible.
255  = -12 : iparam[0] must be equal to 0 or 1.
256  = -13 : nev and which = 'BE' are incompatible.
257  = -9999: Could not build an Arnoldi factorization. iparam[5]
258  returns the size of the current Arnoldi factorization.
259  The user is advised to check that enough workspace
260  and array storage has been allocated.
261 
262  Remarks:
263  1. The converged Ritz values are always returned in ascending
264  algebraic order. The computed Ritz values are approximate
265  eigenvalues of OP. The selection of "which" should be made
266  with this in mind when Mode = 3,4,5. After convergence,
267  approximate eigenvalues of the original problem may be obtained
268  with the ARPACK subroutine seupp.
269  2. If the Ritz vectors corresponding to the converged Ritz values are
270  needed, the user must call seupp immediately following completion
271  of saupp. This is new starting with version 2.1 of ARPACK.
272  3. If M can be factored into a Cholesky factorization M = LL'
273  then Mode = 2 should not be selected. Instead one should use
274  Mode = 1 with OP = inv(L)*A*inv(L'). Appropriate triangular
275  linear systems should be solved with L and L' rather
276  than computing inverses. After convergence, an approximate
277  eigenvector z of the original problem is recovered by solving
278  L'z = x where x is a Ritz vector of OP.
279  4. At present there is no a-priori analysis to guide the selection
280  of ncv relative to nev. The only formal requrement is that
281  ncv > nev. However, it is recommended that ncv >= 2*nev. If many
282  problems of the same type are to be solved, one should experiment
283  with increasing ncv while keeping nev fixed for a given test
284  problem. This will usually decrease the required number of OP*x
285  operations but it also increases the work and storage required to
286  maintain the orthogonal basis vectors. The optimal "cross-over"
287  with respect to CPU time is problem dependent and must be
288  determined empirically.
289  5. If iparam[7] = 2 then in the Reverse commuication interface the
290  user must do the following. When ido = 1, Y = OP * X is to be
291  computed. When iparam[7] = 2 OP = inv(B)*A. After computing A*X
292  the user must overwrite X with A*X. Y is then the solution to the
293  linear set of equations B*Y = A*X.
294  6. When iparam[0] = 0, and ido = 3, the user needs to provide the
295  NP = iparam[8] shifts in locations:
296  1 workl[ipntr[11]]
297  2 workl[ipntr[11]+1]
298  .
299  .
300  .
301  NP workl[ipntr[11]+NP-1].
302  The eigenvalues of the current tridiagonal matrix are located in
303  workl[ipntr[6]] through workl[ipntr[6]+ncv]. They are in the
304  order defined by which. The associated Ritz estimates are located in
305  workl[ipntr[8]], workl[ipntr[8]+1], ... , workl[ipntr[8]+ncv-1].
306 */
307 
308 {
309 
310  F77_FUNC(dsaupd,DSAUPD)(&ido, &bmat, &n, which, &nev, &tol, resid, &ncv,
311  &V[0], &ldv, &iparam[0], &ipntr[0], &workd[0], &workl[0],
312  &lworkl, &info);
313 
314 } // saupp (double).
315 
316 #endif // SAUPP_H