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