tensor-0.1.0
 All Data Structures Namespaces Functions Variables Typedefs Enumerations Enumerator Groups Pages
neupp.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 neupp.h.
24  Interface to ARPACK subroutines dneupd and sneupd.
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 NEUPP_H
36 #define NEUPP_H
37 
38 #include "arpackf.h"
39 
40 inline void neupp(bool rvec, char HowMny, double dr[],
41  double di[], double Z[], int ldz, double sigmar,
42  double sigmai, double workv[], char bmat, int n,
43  char* which, int nev, double tol, double resid[],
44  int ncv, double V[], int ldv, int iparam[],
45  int ipntr[], double workd[], double workl[],
46  int lworkl, int& info)
47 
48 /*
49  c++ version of ARPACK routine dneupd.
50  This subroutine returns the converged approximations to eigenvalues
51  of A*z = lambda*B*z and (optionally):
52 
53  (1) the corresponding approximate eigenvectors,
54  (2) an orthonormal basis for the associated approximate
55  invariant subspace,
56 
57  There is negligible additional cost to obtain eigenvectors. An
58  orthonormal basis is always computed. There is an additional storage cost
59  of n*nev if both are requested (in this case a separate array Z must be
60  supplied).
61  The approximate eigenvalues and eigenvectors of A*z = lambda*B*z
62  are derived from approximate eigenvalues and eigenvectors of
63  of the linear operator OP prescribed by the MODE selection in the
64  call to naupp. naupp must be called before this routine is called.
65  These approximate eigenvalues and vectors are commonly called Ritz
66  values and Ritz vectors respectively. They are referred to as such
67  in the comments that follow. The computed orthonormal basis for the
68  invariant subspace corresponding to these Ritz values is referred to
69  as a Schur basis.
70  See documentation in the header of the subroutine naupp for
71  definition of OP as well as other terms and the relation of computed
72  Ritz values and Ritz vectors of OP with respect to the given problem
73  A*z = lambda*B*z. For a brief description, see definitions of
74  iparam[7], MODE and which in the documentation of naupp.
75 
76  Parameters:
77 
78  rvec (Input) Specifies whether Ritz vectors corresponding to the
79  Ritz value approximations to the eigenproblem A*z = lambda*B*z
80  are computed.
81  rvec = false: Compute Ritz values only.
82  rvec = true : Compute the Ritz vectors or Schur vectors.
83  See Remarks below.
84  HowMny (Input) Specifies the form of the basis for the invariant
85  subspace corresponding to the converged Ritz values that
86  is to be computed.
87  = 'A': Compute nev Ritz vectors;
88  = 'P': Compute nev Schur vectors;
89  dr (Output) Array of dimension nev+1.
90  If iparam[7] = 1,2 or 3 and sigmai=0.0 then on exit: dr
91  contains the real part of the Ritz approximations to the
92  eigenvalues of A*z = lambda*B*z.
93  If iparam[7] = 3, 4 and sigmai is not equal to zero, then on
94  exit: dr contains the real part of the Ritz values of OP
95  computed by naupp. A further computation must be performed by
96  the user to transform the Ritz values computed for OP by naupp
97  to those of the original system A*z = lambda*B*z. See remark 3.
98  di (Output) Array of dimension nev+1.
99  On exit, di contains the imaginary part of the Ritz value
100  approximations to the eigenvalues of A*z = lambda*B*z
101  associated with dr.
102  NOTE: When Ritz values are complex, they will come in complex
103  conjugate pairs. If eigenvectors are requested, the
104  corresponding Ritz vectors will also come in conjugate
105  pairs and the real and imaginary parts of these are
106  represented in two consecutive columns of the array Z
107  (see below).
108  Z (Output) Array of dimension nev*n if rvec = TRUE and HowMny =
109  'A'. if rvec = TRUE. and HowMny = 'A', then the contains
110  approximate eigenvectors (Ritz vectors) corresponding to the
111  NCONV=iparam[5] Ritz values for eigensystem A*z = lambda*B*z.
112  The complex Ritz vector associated with the Ritz value
113  with positive imaginary part is stored in two consecutive
114  columns. The first column holds the real part of the Ritz
115  vector and the second column holds the imaginary part. The
116  Ritz vector associated with the Ritz value with negative
117  imaginary part is simply the complex conjugate of the Ritz
118  vector associated with the positive imaginary part.
119  If rvec = .FALSE. or HowMny = 'P', then Z is not referenced.
120  NOTE: If if rvec = .TRUE. and a Schur basis is not required,
121  the array Z may be set equal to first nev+1 columns of
122  the Arnoldi basis array V computed by naupp. In this
123  case the Arnoldi basis will be destroyed and overwritten
124  with the eigenvector basis.
125  ldz (Input) Dimension of the vectors contained in Z. This
126  parameter MUST be set to n.
127  sigmar (Input) If iparam[7] = 3 or 4, represents the real part of
128  the shift. Not referenced if iparam[7] = 1 or 2.
129  sigmai (Input) If iparam[7] = 3 or 4, represents the imaginary part
130  of the shift. Not referenced if iparam[7] = 1 or 2. See
131  remark 3 below.
132  workv (Workspace) Array of dimension 3*ncv.
133  V (Input/Output) Array of dimension n*ncv+1.
134  Upon Input: V contains the ncv vectors of the Arnoldi basis
135  for OP as constructed by naupp.
136  Upon Output: If rvec = TRUE the first NCONV=iparam[5] columns
137  contain approximate Schur vectors that span the
138  desired invariant subspace. See Remark 2 below.
139  NOTE: If the array Z has been set equal to first nev+1 columns
140  of the array V and rvec = TRUE. and HowMny = 'A', then
141  the Arnoldi basis held by V has been overwritten by the
142  desired Ritz vectors. If a separate array Z has been
143  passed then the first NCONV=iparam[5] columns of V will
144  contain approximate Schur vectors that span the desired
145  invariant subspace.
146  workl (Input / Output) Array of length lworkl+1.
147  workl[1:ncv*ncv+3*ncv] contains information obtained in
148  naupp. They are not changed by neupp.
149  workl[ncv*ncv+3*ncv+1:3*ncv*ncv+6*ncv] holds the real and
150  imaginary part of the untransformed Ritz values, the upper
151  quasi-triangular matrix for H, and the associated matrix
152  representation of the invariant subspace for H.
153  ipntr (Input / Output) Array of length 14. Pointer to mark the
154  starting locations in the workl array for matrices/vectors
155  used by naupp and neupp.
156  ipntr[9]: pointer to the real part of the ncv RITZ values
157  of the original system.
158  ipntr[10]: pointer to the imaginary part of the ncv RITZ
159  values of the original system.
160  ipntr[11]: pointer to the ncv corresponding error bounds.
161  ipntr[12]: pointer to the ncv by ncv upper quasi-triangular
162  Schur matrix for H.
163  ipntr[13]: pointer to the ncv by ncv matrix of eigenvectors
164  of the upper Hessenberg matrix H. Only referenced
165  by neupp if rvec = TRUE. See Remark 2 below.
166  info (Output) Error flag.
167  = 0 : Normal exit.
168  = 1 : The Schur form computed by LAPACK routine dlahqr
169  could not be reordered by LAPACK routine dtrsen.
170  Re-enter subroutine neupp with iparam[5] = ncv and
171  increase the size of the arrays DR and DI to have
172  dimension at least dimension ncv and allocate at least
173  ncv columns for Z. NOTE: Not necessary if Z and V share
174  the same space. Please notify the authors if this error
175  occurs.
176  = -1 : n must be positive.
177  = -2 : nev must be positive.
178  = -3 : ncv must satisfy nev+2 <= ncv <= n.
179  = -5 : which must be one of 'LM','SM','LR','SR','LI','SI'.
180  = -6 : bmat must be one of 'I' or 'G'.
181  = -7 : Length of private work workl array is not sufficient.
182  = -8 : Error return from calculation of a real Schur form.
183  Informational error from LAPACK routine dlahqr.
184  = -9 : Error return from calculation of eigenvectors.
185  Informational error from LAPACK routine dtrevc.
186  = -10: iparam[7] must be 1,2,3,4.
187  = -11: iparam[7] = 1 and bmat = 'G' are incompatible.
188  = -12: HowMny = 'S' not yet implemented
189  = -13: HowMny must be one of 'A' or 'P' if rvec = TRUE.
190  = -14: naupp did not find any eigenvalues to sufficient
191  accuracy.
192 
193  NOTE: The following arguments
194 
195  bmat, n, which, nev, tol, resid, ncv, V, ldv, iparam,
196  ipntr, workd, workl, lworkl, info
197 
198  must be passed directly to neupp following the last call
199  to naupp. These arguments MUST NOT BE MODIFIED between
200  the the last call to naupp and the call to neupp.
201 
202  Remarks
203  1. Currently only HowMny = 'A' and 'P' are implemented.
204  2. Schur vectors are an orthogonal representation for the basis of
205  Ritz vectors. Thus, their numerical properties are often superior.
206  Let X' denote the transpose of X. If rvec = .TRUE. then the
207  relationship A * V[:,1:iparam[5]] = V[:,1:iparam[5]] * T, and
208  V[:,1:iparam[5]]' * V[:,1:iparam[5]] = I are approximately satisfied.
209  Here T is the leading submatrix of order iparam[5] of the real
210  upper quasi-triangular matrix stored workl[ipntr[12]]. That is,
211  T is block upper triangular with 1-by-1 and 2-by-2 diagonal blocks;
212  each 2-by-2 diagonal block has its diagonal elements equal and its
213  off-diagonal elements of opposite sign. Corresponding to each
214  2-by-2 diagonal block is a complex conjugate pair of Ritz values.
215  The real Ritz values are stored on the diagonal of T.
216  3. If iparam[7] = 3 or 4 and sigmai is not equal zero, then the user
217  must form the iparam[5] Rayleigh quotients in order to transform the
218  Ritz values computed by naupp for OP to those of A*z = lambda*B*z.
219  Set rvec = TRUE. and HowMny = 'A', and compute
220  Z[:,I]' * A * Z[:,I] if di[I] = 0.
221  If di[I] is not equal to zero and di[I+1] = - D[I],
222  then the desired real and imaginary parts of the Ritz value are
223  Z[:,I]' * A * Z[:,I] + Z[:,I+1]' * A * Z[:,I+1],
224  Z[:,I]' * A * Z[:,I+1] - Z[:,I+1]' * A * Z[:,I], respectively.
225  Another possibility is to set rvec = .true. and HowMny = 'P' and
226  compute V[:,1:iparam[5]]' * A * V[:,1:iparam[5]] and then an upper
227  quasi-triangular matrix of order iparam[5] is computed. See remark
228  2 above.
229 */
230 
231 {
232 
233  int irvec;
234  logical* iselect;
235  double* iZ;
236 
237  irvec = (int) rvec;
238  iselect = new_atomic logical[ncv];
239  iZ = (Z == NULL) ? &V[0] : Z;
240 
241  F77NAME(dneupd)(&irvec, &HowMny, iselect, dr, di, iZ, &ldz, &sigmar,
242  &sigmai, &workv[0], &bmat, &n, which, &nev, &tol,
243  resid, &ncv, &V[0], &ldv, &iparam[0], &ipntr[0],
244  &workd[0], &workl[0], &lworkl, &info);
245 
246  delete[] iselect;
247 
248 } // neupp (double).
249 
250 #endif // NEUPP_H
251 
252 // Local variables:
253 // mode: c++
254 // fill-column: 80
255 // c-basic-offset: 4