tensor-0.1.0
 All Data Structures Namespaces Functions Variables Typedefs Enumerations Enumerator Groups Pages
ceupp.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 ceupp.h.
24  Interface to ARPACK subroutines zneupd and cneupd.
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 CEUPP_H
36 #define CEUPP_H
37 
38 #include "arpackf.h"
39 
40 inline void ceupp(integer rvec, char HowMny, tensor::cdouble d[],
41  tensor::cdouble Z[], integer ldz, tensor::cdouble sigma,
42  tensor::cdouble workev[], char bmat, integer n, const char* which,
43  integer nev, double tol, tensor::cdouble resid[], integer ncv,
44  tensor::cdouble V[], integer ldv, integer iparam[], integer ipntr[],
45  tensor::cdouble workd[], tensor::cdouble workl[],
46  integer lworkl, double rwork[], integer& info)
47 
48 /*
49  c++ version of ARPACK routine zneupd.
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 caupp. caupp 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 caupp 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 caupp.
75 
76  Parameters:
77 
78  rvec (Input) Specifies whether a basis for the invariant subspace
79  corresponding to the converged Ritz value approximations for
80  the eigenproblem A*z = lambda*B*z is 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  d (Output) Array of dimension nev+1. D contains the Ritz
90  approximations to the eigenvalues lambda for A*z = lambda*B*z.
91  Z (Output) Array of dimension nev*n. If rvec = TRUE. and
92  HowMny = 'A', then Z contains approximate eigenvectors (Ritz
93  vectors) corresponding to the NCONV=iparam[5] Ritz values for
94  eigensystem A*z = lambda*B*z.
95  If rvec = .FALSE. or HowMny = 'P', then Z is not referenced.
96  NOTE: If if rvec = .TRUE. and a Schur basis is not required,
97  the array Z may be set equal to first nev+1 columns of
98  the Arnoldi basis array V computed by caupp. In this
99  case the Arnoldi basis will be destroyed and overwritten
100  with the eigenvector basis.
101  ldz (Input) Dimension of the vectors contained in Z. This
102  parameter MUST be set to n.
103  sigma (Input) If iparam[7] = 3, sigma represents the shift. Not
104  referenced if iparam[7] = 1 or 2.
105  workv (Workspace) Array of dimension 2*ncv.
106  V (Input/Output) Array of dimension n*ncv+1.
107  Upon Input: V contains the ncv vectors of the Arnoldi basis
108  for OP as constructed by caupp.
109  Upon Output: If rvec = TRUE the first NCONV=iparam[5] columns
110  contain approximate Schur vectors that span the
111  desired invariant subspace.
112  NOTE: If the array Z has been set equal to first nev+1 columns
113  of the array V and rvec = TRUE. and HowMny = 'A', then
114  the Arnoldi basis held by V has been overwritten by the
115  desired Ritz vectors. If a separate array Z has been
116  passed then the first NCONV=iparam[5] columns of V will
117  contain approximate Schur vectors that span the desired
118  invariant subspace.
119  workl (Input / Output) Array of length lworkl+1.
120  workl[1:ncv*ncv+3*ncv] contains information obtained in
121  caupp. They are not changed by ceupp.
122  workl[ncv*ncv+3*ncv+1:3*ncv*ncv+4*ncv] holds the untransformed
123  Ritz values, the untransformed error estimates of the Ritz
124  values, the upper triangular matrix for H, and the associated
125  matrix representation of the invariant subspace for H.
126  ipntr (Input / Output) Array of length 14. Pointer to mark the
127  starting locations in the workl array for matrices/vectors
128  used by caupp and ceupp.
129  ipntr[9]: pointer to the ncv RITZ values of the original
130  system.
131  ipntr[11]: pointer to the ncv corresponding error estimates.
132  ipntr[12]: pointer to the ncv by ncv upper triangular
133  Schur matrix for H.
134  ipntr[13]: pointer to the ncv by ncv matrix of eigenvectors
135  of the upper Hessenberg matrix H. Only referenced
136  by ceupp if rvec = TRUE. See Remark 2 below.
137  info (Output) Error flag.
138  = 0 : Normal exit.
139  = 1 : The Schur form computed by LAPACK routine csheqr
140  could not be reordered by LAPACK routine ztrsen.
141  Re-enter subroutine ceupp with iparam[5] = ncv and
142  increase the size of the array D to have
143  dimension at least dimension ncv and allocate at least
144  ncv columns for Z. NOTE: Not necessary if Z and V share
145  the same space. Please notify the authors if this error
146  occurs.
147  = -1 : n must be positive.
148  = -2 : nev must be positive.
149  = -3 : ncv must satisfy nev+1 <= ncv <= n.
150  = -5 : which must be one of 'LM','SM','LR','SR','LI','SI'.
151  = -6 : bmat must be one of 'I' or 'G'.
152  = -7 : Length of private work workl array is not sufficient.
153  = -8 : Error return from LAPACK eigenvalue calculation.
154  This should never happened.
155  = -9 : Error return from calculation of eigenvectors.
156  Informational error from LAPACK routine ztrevc.
157  = -10: iparam[7] must be 1, 2 or 3.
158  = -11: iparam[7] = 1 and bmat = 'G' are incompatible.
159  = -12: HowMny = 'S' not yet implemented.
160  = -13: HowMny must be one of 'A' or 'P' if rvec = TRUE.
161  = -14: caupp did not find any eigenvalues to sufficient
162  accuracy.
163 
164  NOTE: The following arguments
165 
166  bmat, n, which, nev, tol, resid, ncv, V, ldv, iparam,
167  ipntr, workd, workl, lworkl, rwork, info
168 
169  must be passed directly to ceupp following the last call
170  to caupp. These arguments MUST NOT BE MODIFIED between
171  the the last call to caupp and the call to ceupp.
172 
173  Remarks
174  1. Currently only HowMny = 'A' and 'P' are implemented.
175  2. Schur vectors are an orthogonal representation for the basis of
176  Ritz vectors. Thus, their numerical properties are often superior.
177  Let X' denote the transpose of X. If rvec = .TRUE. then the
178  relationship A * V[:,1:iparam[5]] = V[:,1:iparam[5]] * T, and
179  V[:,1:iparam[5]]' * V[:,1:iparam[5]] = I are approximately satisfied.
180  Here T is the leading submatrix of order iparam[5] of the real
181  upper quasi-triangular matrix stored workl[ipntr[12]].
182 */
183 
184 {
185  logical* iselect = new logical[ncv];
186  blas::cdouble* iZ = reinterpret_cast<blas::cdouble*>((Z == NULL) ? V : Z);
187 
188  F77_FUNC(zneupd,ZNEUPD)(&rvec, &HowMny, iselect,
189  reinterpret_cast<blas::cdouble*>(d), iZ, &ldz,
190  reinterpret_cast<blas::cdouble*>(&sigma),
191  reinterpret_cast<blas::cdouble*>(workev),
192  &bmat, &n, which, &nev, &tol,
193  reinterpret_cast<blas::cdouble*>(resid),
194  &ncv, reinterpret_cast<blas::cdouble*>(V),
195  &ldv, &iparam[0], &ipntr[0],
196  reinterpret_cast<blas::cdouble*>(workd),
197  reinterpret_cast<blas::cdouble*>(workl),
198  &lworkl, &rwork[0], &info);
199 
200  delete[] iselect;
201 
202 } // ceupp (cdouble).
203 
204 #endif // CEUPP_H
205 // Local variables:
206 // mode: c++
207 // fill-column: 80
208 // c-basic-offset: 4