tensor-0.1.0
 All Data Structures Namespaces Functions Variables Typedefs Enumerations Enumerator Groups Pages
caupp.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 caupp.h.
24  Interface to ARPACK subroutines znaupd and cnaupd.
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 CAUPP_H
36 #define CAUPP_H
37 
38 #include "arpackf.h"
39 
40 inline void caupp(integer& ido, char bmat, integer n, const char* which, integer nev,
41  double& tol, tensor::cdouble resid[], integer ncv,
42  tensor::cdouble V[], integer ldv, integer iparam[], integer ipntr[],
43  tensor::cdouble workd[], tensor::cdouble workl[],
44  integer lworkl, double rwork[], integer& info)
45 
46 /*
47  c++ version of ARPACK routine znaupd that implements the
48  Reverse communication interface for the Implicitly Restarted Arnoldi
49  iteration. This is intended to be used to find a few eigenpairs of a
50  complex linear operator OP with respect to a semi-inner product defined
51  by a hermitian positive semi-definite real matrix B. B may be the
52  identity matrix. NOTE: if both OP and B are real, then naupp should
53  be used.
54 
55  The computed approximate eigenvalues are called Ritz values and
56  the corresponding approximate eigenvectors are called Ritz vectors.
57 
58  caupp is usually called iteratively to solve one of the
59  following problems:
60 
61  Mode 1: A*x = lambda*x.
62  ===> OP = A and B = I.
63 
64  Mode 2: A*x = lambda*M*x, M symmetric positive definite
65  ===> OP = inv[M]*A and B = M.
66  ===> (If M can be factored see remark 3 below)
67 
68  Mode 3: A*x = lambda*M*x, M symmetric semi-definite
69  ===> OP = inv[A - sigma*M]*M and B = M.
70  ===> shift-and-invert mode
71  If OP*x = amu*x, then lambda = sigma + 1/amu.
72 
73 
74  NOTE: The action of w <- inv[A - sigma*M]*v or w <- inv[M]*v
75  should be accomplished either by a direct method
76  using a sparse matrix factorization and solving
77 
78  [A - sigma*M]*w = v or M*w = v,
79 
80  or through an iterative method for solving these systems. If
81  an iterative method is used, the convergence test must be more
82  stringent than the accuracy requirements for the eigenvalue
83  approximations.
84 
85  Parameters:
86 
87  ido (Input / Output) Reverse communication flag. ido must be
88  zero on the first call to caupp. ido will be set
89  internally to indicate the type of operation to be
90  performed. Control is then given back to the calling
91  routine which has the responsibility to carry out the
92  requested operation and call caupp with the result. The
93  operand is given in workd[ipntr[0]], the result must be
94  put in workd[ipntr[2]].
95  ido = 0: first call to the reverse communication interface.
96  ido = -1: compute Y = OP * X where
97  ipntr[0] is the pointer into workd for X,
98  ipntr[2] is the pointer into workd for Y.
99  This is for the initialization phase to force the
100  starting vector into the range of OP.
101  ido = 1: compute Y = OP * X where
102  ipntr[0] is the pointer into workd for X,
103  ipntr[2] is the pointer into workd for Y.
104  In mode 3 and 4, the vector B * X is already
105  available in workd[ipntr[3]]. It does not
106  need to be recomputed in forming OP * X.
107  ido = 2: compute Y = B * X where
108  ipntr[0] is the pointer into workd for X,
109  ipntr[2] is the pointer into workd for Y.
110  ido = 3: compute the iparam[8] real and imaginary parts
111  of the shifts where inptr[14] is the pointer
112  into workl for placing the shifts. See Remark
113  5 below.
114  ido = 99: done.
115  bmat (Input) bmat specifies the type of the matrix B that defines
116  the semi-inner product for the operator OP.
117  bmat = 'I' -> standard eigenvalue problem A*x = lambda*x;
118  bmat = 'G' -> generalized eigenvalue problem A*x = lambda*M*x.
119  n (Input) Dimension of the eigenproblem.
120  nev (Input) Number of eigenvalues to be computed. 0 < nev <= n-1.
121  which (Input) Specify which of the Ritz values of OP to compute.
122  'LM' - compute the nev eigenvalues of largest magnitude.
123  'SM' - compute the nev eigenvalues of smallest magnitude.
124  'LR' - compute the nev eigenvalues of largest real part.
125  'SR' - compute the nev eigenvalues of smallest real part.
126  'LI' - compute the nev eigenvalues of largest imaginary part.
127  'SI' - compute the nev eigenvalues of smallest imaginary part.
128  tol (Input) Stopping criterion: the relative accuracy of the
129  Ritz value is considered acceptable if BOUNDS[i] <=
130  tol*abs(RITZ[i]),where ABS(RITZ[i]) is the magnitude when
131  RITZ[i] is complex. If tol<=0.0 is passed, the machine
132  precision as computed by the LAPACK auxiliary subroutine
133  _LAMCH is used.
134  resid (Input / Output) Array of length n.
135  On input:
136  If info==0, a random initial residual vector is used.
137  If info!=0, resid contains the initial residual vector,
138  possibly from a previous run.
139  On output:
140  resid contains the final residual vector.
141  ncv (Input) Number of Arnoldi vectors that are generated at each
142  iteration. After the startup phase in which nev Arnoldi
143  vectors are generated, the algorithm generates ncv-nev
144  Arnoldi vectors at each subsequent update iteration. Most of
145  the cost in generating each Arnoldi vector is in the
146  matrix-vector product OP*x.
147  NOTE: ncv must satisfy nev+1 <= ncv <= n.
148  V (Output) Array of length ncv*n+1. V contains the ncv Arnoldi
149  basis vectors. The first element V[0] is never referenced.
150  ldv (Input) Dimension of the basis vectors contained in V. This
151  parameter MUST be set to n.
152  iparam (Input / Output) Array of length 12.
153  iparam[0] = ISHIFT: method for selecting the implicit shifts.
154  The shifts selected at each iteration are used to restart
155  the Arnoldi iteration in an implicit fashion.
156  -------------------------------------------------------------
157  ISHIFT = 0: the shifts are to be provided by the user via
158  reverse communication. The ncv eigenvalues of
159  the Hessenberg matrix H are returned in the part
160  of workl array corresponding to RITZ.
161  ISHIFT = 1: exact shifts with respect to the current
162  Hessenberg matrix H. This is equivalent to
163  restarting the iteration from the beginning
164  after updating the starting vector with a linear
165  combination of Ritz vectors associated with the
166  "wanted" eigenvalues.
167  ISHIFT = 2: other choice of internal shift to be defined.
168  -------------------------------------------------------------
169  iparam[2] is no longer referenced.
170  iparam[3] = MXITER
171  On INPUT: maximum number of Arnoldi update iterations allowed.
172  On OUTPUT: actual number of Arnoldi update iterations taken.
173  iparam[4] = NB: blocksize to be used in the recurrence.
174  The code currently works only for NB = 1.
175  iparam[5] = NCONV: number of "converged" Ritz values.
176  This represents the number of Ritz values that satisfy
177  the convergence criterion.
178  iparam[6] is no longer referenced.
179  iparam[7] = MODE. On input determines what type of
180  eigenproblem is being solved. Must be 1, 2 or 3.
181  iparam[8] = NP. When ido = 3 and the user provides shifts
182  through reverse communication (iparam[0]=0), caupp returns
183  NP, the number of shifts the user is to provide.
184  0 < NP <=ncv-nev. See Remark 5 below.
185  iparam[9] = total number of OP*x operations.
186  iparam[10] = total number of B*x operations if bmat='G'.
187  iparam[11] = total number of steps of re-orthogonalization.
188  ipntr (Output) Array of length 15. Pointer to mark the starting
189  locations in the workd and workl arrays for matrices/vectors
190  used by the Arnoldi iteration.
191  ipntr[0] : pointer to the current operand vector X in workd.
192  ipntr[2] : pointer to the current result vector Y in workd.
193  ipntr[3] : pointer to the vector B * X in workd when used in
194  the shift-and-invert mode.
195  ipntr[4] : pointer to the next available location in workl
196  that is untouched by the program.
197  ipntr[5] : pointer to the ncv by ncv upper Hessenberg matrix
198  H in workl.
199  ipntr[6] : pointer to the ritz value array RITZ.
200  ipntr[7] : pointer to the (projected) ritz vector array Q.
201  ipntr[8] : pointer to the error BOUNDS array in workl.
202  ipntr[14]: pointer to the NP shifts in workl. See Remark 5.
203  Note: ipntr[9:13] is only referenced by ceupp. See Remark 2.
204  ipntr[9] : pointer to the ncv RITZ values of the
205  original system.
206  ipntr[10]: Not Used
207  ipntr[11]: pointer to the ncv corresponding error bounds.
208  ipntr[12]: pointer to the ncv by ncv upper triangular
209  Schur matrix for H.
210  ipntr[13]: pointer to the ncv by ncv matrix of eigenvectors
211  of the upper Hessenberg matrix H. Only referenced by
212  ceupp if RVEC = true. See Remark 2 below.
213  workd (Input / Output) Array of length 3*n+1.
214  Distributed array to be used in the basic Arnoldi iteration
215  for reverse communication. The user should not use workd as
216  temporary workspace during the iteration.
217  workl (Output) Array of length lworkl+1. Private (replicated) array
218  on each PE or array allocated on the front end.
219  lworkl (Input) lworkl must be at least 3*ncv*ncv+5*ncv.
220  RWORK (Workspace) Array of length ncv. Private (replicated) array on
221  each PE or array allocated on the front end.
222  info (Input / Output) On input, if info = 0, a randomly initial
223  residual vector is used, otherwise resid contains the initial
224  residual vector, possibly from a previous run.
225  On output, info works as a error flag:
226  = 0 : Normal exit.
227  = 1 : Maximum number of iterations taken. All possible
228  eigenvalues of OP has been found. iparam[5]
229  returns the number of wanted converged Ritz values.
230  = 3 : No shifts could be applied during a cycle of the
231  Implicitly restarted Arnoldi iteration. One
232  possibility is to increase the size of ncv relative
233  to nev. See remark 4 below.
234  = -1 : n must be positive.
235  = -2 : nev must be positive.
236  = -3 : ncv must satisfy nev+1 <= ncv <= n.
237  = -4 : The maximum number of Arnoldi update iterations
238  allowed must be greater than zero.
239  = -5 : which must be one of 'LM','SM','LR','SR','LI','SI'.
240  = -6 : bmat must be one of 'I' or 'G'.
241  = -7 : Length of private work array is not sufficient.
242  = -8 : Error return from LAPACK eigenvalue calculation.
243  = -9 : Starting vector is zero.
244  = -10 : iparam[7] must be 1, 2 or 3.
245  = -11 : iparam[7] = 1 and bmat = 'G' are incompatible.
246  = -12 : iparam[0] must be equal to 0 or 1.
247  = -13 : nev and which = 'BE' are incompatible.
248  = -9999: Could not build an Arnoldi factorization. iparam[5]
249  returns the size of the current Arnoldi factorization.
250  The user is advised to check that enough workspace
251  and array storage has been allocated.
252 
253  Remarks:
254  1. The computed Ritz values are approximate eigenvalues of OP. The
255  selection of "which" should be made with this in mind when using
256  Mode = 3. When operating in Mode = 3 setting which = 'LM' will
257  compute the nev eigenvalues of the original problem that are
258  closest to the shift sigma . After convergence, approximate
259  eigenvalues of the original problem may be obtained with the
260  ARPACK subroutine ceupp.
261  2. If a basis for the invariant subspace corresponding to the converged
262  Ritz values is needed, the user must call ceupp immediately following
263  completion of caupp. This is new starting with release 2 of ARPACK.
264  3. If M can be factored into a Cholesky factorization M = LL'
265  then Mode = 2 should not be selected. Instead one should use
266  Mode = 1 with OP = inv(L)*A*inv(L'). Appropriate triangular
267  linear systems should be solved with L and L' rather
268  than computing inverses. After convergence, an approximate
269  eigenvector z of the original problem is recovered by solving
270  L'z = x where x is a Ritz vector of OP.
271  4. At present there is no a-priori analysis to guide the selection
272  of ncv relative to nev. The only formal requrement is that ncv
273  >= nev + 1. However, it is recommended that ncv >= 2*nev. If many
274  problems of the same type are to be solved, one should experiment
275  with increasing ncv while keeping nev fixed for a given test
276  problem. This will usually decrease the required number of OP*x
277  operations but it also increases the work and storage required to
278  maintain the orthogonal basis vectors. The optimal "cross-over"
279  with respect to CPU time is problem dependent and must be
280  determined empirically.
281  5. When iparam[0] = 0, and ido = 3, the user needs to provide the
282  NP = iparam[8] complex shifts in locations
283  workl[ipntr[14]], workl[ipntr[14]+1], ... , workl[ipntr[14]+NP].
284  Eigenvalues of the current upper Hessenberg matrix are located in
285  workl[ipntr[6]] through workl[ipntr[6]+ncv-1]. They are ordered
286  according to the order defined by "which". The associated Ritz
287  estimates are located in workl[ipntr[8]], workl[ipntr[8]+1], ...,
288  workl[ipntr[8]+ncv-1].
289 
290  References:
291  1. D.C. Sorensen, "Implicit Application of Polynomial Filters in
292  a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992),
293  pp 357-385.
294  2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly
295  Restarted Arnoldi Iteration", Rice University Technical Report
296  TR95-13, Department of Computational and Applied Mathematics.
297  3. B.N. Parlett & Y. Saad, "_Complex_ Shift and Invert Strategies for
298  Double precision Matrices", Linear Algebra and its Applications,
299  vol 88/89, pp 575-595, (1987).
300 */
301 
302 {
303  F77_FUNC(znaupd,ZNAUPD)(&ido, &bmat, &n, which, &nev, &tol,
304  reinterpret_cast<blas::cdouble*>(resid), &ncv,
305  reinterpret_cast<blas::cdouble*>(V), &ldv,
306  iparam, ipntr,
307  reinterpret_cast<blas::cdouble*>(workd),
308  reinterpret_cast<blas::cdouble*>(workl),
309  &lworkl, rwork, &info);
310 } // caupp (cdouble).
311 
312 #endif // CAUPP_H