1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/libraries/libmbwrap/arptest.cc,v 1.12 2017/01/12 14:44:25 masarati Exp $ */
2 /*
3 * MBDyn (C) is a multibody analysis code.
4 * http://www.mbdyn.org
5 *
6 * Copyright (C) 1996-2017
7 *
8 * Pierangelo Masarati <masarati@aero.polimi.it>
9 * Paolo Mantegazza <mantegazza@aero.polimi.it>
10 *
11 * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
12 * via La Masa, 34 - 20156 Milano, Italy
13 * http://www.aero.polimi.it
14 *
15 * Changing this copyright notice is forbidden.
16 *
17 * This program is free software; you can redistribute it and/or modify
18 * it under the terms of the GNU General Public License as published by
19 * the Free Software Foundation (version 2 of the License).
20 *
21 *
22 * This program is distributed in the hope that it will be useful,
23 * but WITHOUT ANY WARRANTY; without even the implied warranty of
24 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 * GNU General Public License for more details.
26 *
27 * You should have received a copy of the GNU General Public License
28 * along with this program; if not, write to the Free Software
29 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
30 */
31
32 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
33
34 #include <iostream>
35 #include <iomanip>
36
37 #include <vector>
38
39 #include "ac/f2c.h"
40 #include "ac/arpack.h"
41
42 #include "solman.h"
43 #include "naivewrap.h"
44
45 /* from dnaupd.f distributed with arpack:
46
47 c\Description:
48 c Reverse communication interface for the Implicitly Restarted Arnoldi
49 c iteration. This subroutine computes approximations to a few eigenpairs
50 c of a linear operator "OP" with respect to a semi-inner product defined by
51 c a symmetric positive semi-definite real matrix B. B may be the identity
52 c matrix. NOTE: If the linear operator "OP" is real and symmetric
53 c with respect to the real positive semi-definite symmetric matrix B,
54 c i.e. B*OP = (OP`)*B, then subroutine dsaupd should be used instead.
55 c
56 c The computed approximate eigenvalues are called Ritz values and
57 c the corresponding approximate eigenvectors are called Ritz vectors.
58 c
59 c dnaupd is usually called iteratively to solve one of the
60 c following problems:
61 c
62 c Mode 1: A*x = lambda*x.
63 c ===> OP = A and B = I.
64 c
65 c Mode 2: A*x = lambda*M*x, M symmetric positive definite
66 c ===> OP = inv[M]*A and B = M.
67 c ===> (If M can be factored see remark 3 below)
68 c
69 c Mode 3: A*x = lambda*M*x, M symmetric semi-definite
70 c ===> OP = Real_Part{ inv[A - sigma*M]*M } and B = M.
71 c ===> shift-and-invert mode (in real arithmetic)
72 c If OP*x = amu*x, then
73 c amu = 1/2 * [ 1/(lambda-sigma) + 1/(lambda-conjg(sigma)) ].
74 c Note: If sigma is real, i.e. imaginary part of sigma is zero;
75 c Real_Part{ inv[A - sigma*M]*M } == inv[A - sigma*M]*M
76 c amu == 1/(lambda-sigma).
77 c
78 c Mode 4: A*x = lambda*M*x, M symmetric semi-definite
79 c ===> OP = Imaginary_Part{ inv[A - sigma*M]*M } and B = M.
80 c ===> shift-and-invert mode (in real arithmetic)
81 c If OP*x = amu*x, then
82 c amu = 1/2i * [ 1/(lambda-sigma) - 1/(lambda-conjg(sigma)) ].
83 c
84 c Both mode 3 and 4 give the same enhancement to eigenvalues close to
85 c the (complex) shift sigma. However, as lambda goes to infinity,
86 c the operator OP in mode 4 dampens the eigenvalues more strongly than
87 c does OP defined in mode 3.
88 c
89 c NOTE: The action of w <- inv[A - sigma*M]*v or w <- inv[M]*v
90 c should be accomplished either by a direct method
91 c using a sparse matrix factorization and solving
92 c
93 c [A - sigma*M]*w = v or M*w = v,
94 c
95 c or through an iterative method for solving these
96 c systems. If an iterative method is used, the
97 c convergence test must be more stringent than
98 c the accuracy requirements for the eigenvalue
99 c approximations.
100 c
101 c\Usage:
102 c call dnaupd
103 c ( IDO, BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM,
104 c IPNTR, WORKD, WORKL, LWORKL, INFO )
105 c
106 c\Arguments
107 c IDO Integer. (INPUT/OUTPUT)
108 c Reverse communication flag. IDO must be zero on the first
109 c call to dnaupd . IDO will be set internally to
110 c indicate the type of operation to be performed. Control is
111 c then given back to the calling routine which has the
112 c responsibility to carry out the requested operation and call
113 c dnaupd with the result. The operand is given in
114 c WORKD(IPNTR(1)), the result must be put in WORKD(IPNTR(2)).
115 c -------------------------------------------------------------
116 c IDO = 0: first call to the reverse communication interface
117 c IDO = -1: compute Y = OP * X where
118 c IPNTR(1) is the pointer into WORKD for X,
119 c IPNTR(2) is the pointer into WORKD for Y.
120 c This is for the initialization phase to force the
121 c starting vector into the range of OP.
122 c IDO = 1: compute Y = OP * X where
123 c IPNTR(1) is the pointer into WORKD for X,
124 c IPNTR(2) is the pointer into WORKD for Y.
125 c In mode 3 and 4, the vector B * X is already
126 c available in WORKD(ipntr(3)). It does not
127 c need to be recomputed in forming OP * X.
128 c IDO = 2: compute Y = B * X where
129 c IPNTR(1) is the pointer into WORKD for X,
130 c IPNTR(2) is the pointer into WORKD for Y.
131 c IDO = 3: compute the IPARAM(8) real and imaginary parts
132 c of the shifts where INPTR(14) is the pointer
133 c into WORKL for placing the shifts. See Remark
134 c 5 below.
135 c IDO = 99: done
136 c -------------------------------------------------------------
137 c
138 c BMAT Character*1. (INPUT)
139 c BMAT specifies the type of the matrix B that defines the
140 c semi-inner product for the operator OP.
141 c BMAT = 'I' -> standard eigenvalue problem A*x = lambda*x
142 c BMAT = 'G' -> generalized eigenvalue problem A*x = lambda*B*x
143 c
144 c N Integer. (INPUT)
145 c Dimension of the eigenproblem.
146 c
147 c WHICH Character*2. (INPUT)
148 c 'LM' -> want the NEV eigenvalues of largest magnitude.
149 c 'SM' -> want the NEV eigenvalues of smallest magnitude.
150 c 'LR' -> want the NEV eigenvalues of largest real part.
151 c 'SR' -> want the NEV eigenvalues of smallest real part.
152 c 'LI' -> want the NEV eigenvalues of largest imaginary part.
153 c 'SI' -> want the NEV eigenvalues of smallest imaginary part.
154 c
155 c NEV Integer. (INPUT)
156 c Number of eigenvalues of OP to be computed. 0 < NEV < N-1.
157 c
158 c TOL Double precision scalar. (INPUT)
159 c Stopping criterion: the relative accuracy of the Ritz value
160 c is considered acceptable if BOUNDS(I) .LE. TOL*ABS(RITZ(I))
161 c where ABS(RITZ(I)) is the magnitude when RITZ(I) is complex.
162 c DEFAULT = DLAMCH ('EPS') (machine precision as computed
163 c by the LAPACK auxiliary subroutine DLAMCH ).
164 c
165 c RESID Double precision array of length N. (INPUT/OUTPUT)
166 c On INPUT:
167 c If INFO .EQ. 0, a random initial residual vector is used.
168 c If INFO .NE. 0, RESID contains the initial residual vector,
169 c possibly from a previous run.
170 c On OUTPUT:
171 c RESID contains the final residual vector.
172 c
173 c NCV Integer. (INPUT)
174 c Number of columns of the matrix V. NCV must satisfy the two
175 c inequalities 2 <= NCV-NEV and NCV <= N.
176 c This will indicate how many Arnoldi vectors are generated
177 c at each iteration. After the startup phase in which NEV
178 c Arnoldi vectors are generated, the algorithm generates
179 c approximately NCV-NEV Arnoldi vectors at each subsequent update
180 c iteration. Most of the cost in generating each Arnoldi vector is
181 c in the matrix-vector operation OP*x.
182 c NOTE: 2 <= NCV-NEV in order that complex conjugate pairs of Ritz
183 c values are kept together. (See remark 4 below)
184 c
185 c V Double precision array N by NCV. (OUTPUT)
186 c Contains the final set of Arnoldi basis vectors.
187 c
188 c LDV Integer. (INPUT)
189 c Leading dimension of V exactly as declared in the calling program.
190 c
191 c IPARAM Integer array of length 11. (INPUT/OUTPUT)
192 c IPARAM(1) = ISHIFT: method for selecting the implicit shifts.
193 c The shifts selected at each iteration are used to restart
194 c the Arnoldi iteration in an implicit fashion.
195 c -------------------------------------------------------------
196 c ISHIFT = 0: the shifts are provided by the user via
197 c reverse communication. The real and imaginary
198 c parts of the NCV eigenvalues of the Hessenberg
199 c matrix H are returned in the part of the WORKL
200 c array corresponding to RITZR and RITZI. See remark
201 c 5 below.
202 c ISHIFT = 1: exact shifts with respect to the current
203 c Hessenberg matrix H. This is equivalent to
204 c restarting the iteration with a starting vector
205 c that is a linear combination of approximate Schur
206 c vectors associated with the "wanted" Ritz values.
207 c -------------------------------------------------------------
208 c
209 c IPARAM(2) = No longer referenced.
210 c
211 c IPARAM(3) = MXITER
212 c On INPUT: maximum number of Arnoldi update iterations allowed.
213 c On OUTPUT: actual number of Arnoldi update iterations taken.
214 c
215 c IPARAM(4) = NB: blocksize to be used in the recurrence.
216 c The code currently works only for NB = 1.
217 c
218 c IPARAM(5) = NCONV: number of "converged" Ritz values.
219 c This represents the number of Ritz values that satisfy
220 c the convergence criterion.
221 c
222 c IPARAM(6) = IUPD
223 c No longer referenced. Implicit restarting is ALWAYS used.
224 c
225 c IPARAM(7) = MODE
226 c On INPUT determines what type of eigenproblem is being solved.
227 c Must be 1,2,3,4; See under \Description of dnaupd for the
228 c four modes available.
229 c
230 c IPARAM(8) = NP
231 c When ido = 3 and the user provides shifts through reverse
232 c communication (IPARAM(1)=0), dnaupd returns NP, the number
233 c of shifts the user is to provide. 0 < NP <=NCV-NEV. See Remark
234 c 5 below.
235 c
236 c IPARAM(9) = NUMOP, IPARAM(10) = NUMOPB, IPARAM(11) = NUMREO,
237 c OUTPUT: NUMOP = total number of OP*x operations,
238 c NUMOPB = total number of B*x operations if BMAT='G',
239 c NUMREO = total number of steps of re-orthogonalization.
240 c
241 c IPNTR Integer array of length 14. (OUTPUT)
242 c Pointer to mark the starting locations in the WORKD and WORKL
243 c arrays for matrices/vectors used by the Arnoldi iteration.
244 c -------------------------------------------------------------
245 c IPNTR(1): pointer to the current operand vector X in WORKD.
246 c IPNTR(2): pointer to the current result vector Y in WORKD.
247 c IPNTR(3): pointer to the vector B * X in WORKD when used in
248 c the shift-and-invert mode.
249 c IPNTR(4): pointer to the next available location in WORKL
250 c that is untouched by the program.
251 c IPNTR(5): pointer to the NCV by NCV upper Hessenberg matrix
252 c H in WORKL.
253 c IPNTR(6): pointer to the real part of the ritz value array
254 c RITZR in WORKL.
255 c IPNTR(7): pointer to the imaginary part of the ritz value array
256 c RITZI in WORKL.
257 c IPNTR(8): pointer to the Ritz estimates in array WORKL associated
258 c with the Ritz values located in RITZR and RITZI in WORKL.
259 c
260 c IPNTR(14): pointer to the NP shifts in WORKL. See Remark 5 below.
261 c
262 c Note: IPNTR(9:13) is only referenced by dneupd . See Remark 2 below.
263 c
264 c IPNTR(9): pointer to the real part of the NCV RITZ values of the
265 c original system.
266 c IPNTR(10): pointer to the imaginary part of the NCV RITZ values of
267 c the original system.
268 c IPNTR(11): pointer to the NCV corresponding error bounds.
269 c IPNTR(12): pointer to the NCV by NCV upper quasi-triangular
270 c Schur matrix for H.
271 c IPNTR(13): pointer to the NCV by NCV matrix of eigenvectors
272 c of the upper Hessenberg matrix H. Only referenced by
273 c dneupd if RVEC = .TRUE. See Remark 2 below.
274 c -------------------------------------------------------------
275 c
276 c WORKD Double precision work array of length 3*N. (REVERSE COMMUNICATION)
277 c Distributed array to be used in the basic Arnoldi iteration
278 c for reverse communication. The user should not use WORKD
279 c as temporary workspace during the iteration. Upon termination
280 c WORKD(1:N) contains B*RESID(1:N). If an invariant subspace
281 c associated with the converged Ritz values is desired, see remark
282 c 2 below, subroutine dneupd uses this output.
283 c See Data Distribution Note below.
284 c
285 c WORKL Double precision work array of length LWORKL. (OUTPUT/WORKSPACE)
286 c Private (replicated) array on each PE or array allocated on
287 c the front end. See Data Distribution Note below.
288 c
289 c LWORKL Integer. (INPUT)
290 c LWORKL must be at least 3*NCV**2 + 6*NCV.
291 c
292 c INFO Integer. (INPUT/OUTPUT)
293 c If INFO .EQ. 0, a randomly initial residual vector is used.
294 c If INFO .NE. 0, RESID contains the initial residual vector,
295 c possibly from a previous run.
296 c Error flag on output.
297 c = 0: Normal exit.
298 c = 1: Maximum number of iterations taken.
299 c All possible eigenvalues of OP has been found. IPARAM(5)
300 c returns the number of wanted converged Ritz values.
301 c = 2: No longer an informational error. Deprecated starting
302 c with release 2 of ARPACK.
303 c = 3: No shifts could be applied during a cycle of the
304 c Implicitly restarted Arnoldi iteration. One possibility
305 c is to increase the size of NCV relative to NEV.
306 c See remark 4 below.
307 c = -1: N must be positive.
308 c = -2: NEV must be positive.
309 c = -3: NCV-NEV >= 2 and less than or equal to N.
310 c = -4: The maximum number of Arnoldi update iteration
311 c must be greater than zero.
312 c = -5: WHICH must be one of 'LM', 'SM', 'LR', 'SR', 'LI', 'SI'
313 c = -6: BMAT must be one of 'I' or 'G'.
314 c = -7: Length of private work array is not sufficient.
315 c = -8: Error return from LAPACK eigenvalue calculation;
316 c = -9: Starting vector is zero.
317 c = -10: IPARAM(7) must be 1,2,3,4.
318 c = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatable.
319 c = -12: IPARAM(1) must be equal to 0 or 1.
320 c = -9999: Could not build an Arnoldi factorization.
321 c IPARAM(5) returns the size of the current Arnoldi
322 c factorization.
323 c
324 c\Remarks
325 c 1. The computed Ritz values are approximate eigenvalues of OP. The
326 c selection of WHICH should be made with this in mind when
327 c Mode = 3 and 4. After convergence, approximate eigenvalues of the
328 c original problem may be obtained with the ARPACK subroutine dneupd .
329 c
330 c 2. If a basis for the invariant subspace corresponding to the converged Ritz
331 c values is needed, the user must call dneupd immediately following
332 c completion of dnaupd . This is new starting with release 2 of ARPACK.
333 c
334 c 3. If M can be factored into a Cholesky factorization M = LL`
335 c then Mode = 2 should not be selected. Instead one should use
336 c Mode = 1 with OP = inv(L)*A*inv(L`). Appropriate triangular
337 c linear systems should be solved with L and L` rather
338 c than computing inverses. After convergence, an approximate
339 c eigenvector z of the original problem is recovered by solving
340 c L`z = x where x is a Ritz vector of OP.
341 c
342 c 4. At present there is no a-priori analysis to guide the selection
343 c of NCV relative to NEV. The only formal requrement is that NCV > NEV + 2.
344 c However, it is recommended that NCV .ge. 2*NEV+1. If many problems of
345 c the same type are to be solved, one should experiment with increasing
346 c NCV while keeping NEV fixed for a given test problem. This will
347 c usually decrease the required number of OP*x operations but it
348 c also increases the work and storage required to maintain the orthogonal
349 c basis vectors. The optimal "cross-over" with respect to CPU time
350 c is problem dependent and must be determined empirically.
351 c See Chapter 8 of Reference 2 for further information.
352 c
353 c 5. When IPARAM(1) = 0, and IDO = 3, the user needs to provide the
354 c NP = IPARAM(8) real and imaginary parts of the shifts in locations
355 c real part imaginary part
356 c ----------------------- --------------
357 c 1 WORKL(IPNTR(14)) WORKL(IPNTR(14)+NP)
358 c 2 WORKL(IPNTR(14)+1) WORKL(IPNTR(14)+NP+1)
359 c . .
360 c . .
361 c . .
362 c NP WORKL(IPNTR(14)+NP-1) WORKL(IPNTR(14)+2*NP-1).
363 c
364 c Only complex conjugate pairs of shifts may be applied and the pairs
365 c must be placed in consecutive locations. The real part of the
366 c eigenvalues of the current upper Hessenberg matrix are located in
367 c WORKL(IPNTR(6)) through WORKL(IPNTR(6)+NCV-1) and the imaginary part
368 c in WORKL(IPNTR(7)) through WORKL(IPNTR(7)+NCV-1). They are ordered
369 c according to the order defined by WHICH. The complex conjugate
370 c pairs are kept together and the associated Ritz estimates are located in
371 c WORKL(IPNTR(8)), WORKL(IPNTR(8)+1), ... , WORKL(IPNTR(8)+NCV-1).
372 c
373 c-----------------------------------------------------------------------
374 c
375 c\Data Distribution Note:
376 c
377 c Fortran-D syntax:
378 c ================
379 c Double precision resid(n), v(ldv,ncv), workd(3*n), workl(lworkl)
380 c decompose d1(n), d2(n,ncv)
381 c align resid(i) with d1(i)
382 c align v(i,j) with d2(i,j)
383 c align workd(i) with d1(i) range (1:n)
384 c align workd(i) with d1(i-n) range (n+1:2*n)
385 c align workd(i) with d1(i-2*n) range (2*n+1:3*n)
386 c distribute d1(block), d2(block,:)
387 c replicated workl(lworkl)
388 c
389 c Cray MPP syntax:
390 c ===============
391 c Double precision resid(n), v(ldv,ncv), workd(n,3), workl(lworkl)
392 c shared resid(block), v(block,:), workd(block,:)
393 c replicated workl(lworkl)
394 c
395 c CM2/CM5 syntax:
396 c ==============
397 c
398 c-----------------------------------------------------------------------
399 c
400 c include 'ex-nonsym.doc'
401 c
402 c-----------------------------------------------------------------------
403 c
404 c\BeginLib
405 c
406 c\Local variables:
407 c xxxxxx real
408 c
409 c\References:
410 c 1. D.C. Sorensen, "Implicit Application of Polynomial Filters in
411 c a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992),
412 c pp 357-385.
413 c 2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly
414 c Restarted Arnoldi Iteration", Rice University Technical Report
415 c TR95-13, Department of Computational and Applied Mathematics.
416 c 3. B.N. Parlett & Y. Saad, "Complex Shift and Invert Strategies for
417 c Real Matrices", Linear Algebra and its Applications, vol 88/89,
418 c pp 575-595, (1987).
419 c
420 c\Routines called:
421 c dnaup2 ARPACK routine that implements the Implicitly Restarted
422 c Arnoldi Iteration.
423 c ivout ARPACK utility routine that prints integers.
424 c second ARPACK utility routine for timing.
425 c dvout ARPACK utility routine that prints vectors.
426 c dlamch LAPACK routine that determines machine constants.
427 c
428 c\Author
429 c Danny Sorensen Phuong Vu
430 c Richard Lehoucq CRPC / Rice University
431 c Dept. of Computational & Houston, Texas
432 c Applied Mathematics
433 c Rice University
434 c Houston, Texas
435 c
436 c\Revision history:
437 c 12/16/93: Version '1.1'
438 c
439 c\SCCS Information: @(#)
440 c FILE: naupd.F SID: 2.8 DATE OF SID: 04/10/01 RELEASE: 2
441 */
442
443 /* from dneupd.f distributed with arpack:
444 c\Name: dneupd
445 c
446 c\Description:
447 c
448 c This subroutine returns the converged approximations to eigenvalues
449 c of A*z = lambda*B*z and (optionally):
450 c
451 c (1) The corresponding approximate eigenvectors;
452 c
453 c (2) An orthonormal basis for the associated approximate
454 c invariant subspace;
455 c
456 c (3) Both.
457 c
458 c There is negligible additional cost to obtain eigenvectors. An orthonormal
459 c basis is always computed. There is an additional storage cost of n*nev
460 c if both are requested (in this case a separate array Z must be supplied).
461 c
462 c The approximate eigenvalues and eigenvectors of A*z = lambda*B*z
463 c are derived from approximate eigenvalues and eigenvectors of
464 c of the linear operator OP prescribed by the MODE selection in the
465 c call to DNAUPD . DNAUPD must be called before this routine is called.
466 c These approximate eigenvalues and vectors are commonly called Ritz
467 c values and Ritz vectors respectively. They are referred to as such
468 c in the comments that follow. The computed orthonormal basis for the
469 c invariant subspace corresponding to these Ritz values is referred to as a
470 c Schur basis.
471 c
472 c See documentation in the header of the subroutine DNAUPD for
473 c definition of OP as well as other terms and the relation of computed
474 c Ritz values and Ritz vectors of OP with respect to the given problem
475 c A*z = lambda*B*z. For a brief description, see definitions of
476 c IPARAM(7), MODE and WHICH in the documentation of DNAUPD .
477 c
478 c\Usage:
479 c call dneupd
480 c ( RVEC, HOWMNY, SELECT, DR, DI, Z, LDZ, SIGMAR, SIGMAI, WORKEV, BMAT,
481 c N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM, IPNTR, WORKD, WORKL,
482 c LWORKL, INFO )
483 c
484 c\Arguments:
485 c RVEC LOGICAL (INPUT)
486 c Specifies whether a basis for the invariant subspace corresponding
487 c to the converged Ritz value approximations for the eigenproblem
488 c A*z = lambda*B*z is computed.
489 c
490 c RVEC = .FALSE. Compute Ritz values only.
491 c
492 c RVEC = .TRUE. Compute the Ritz vectors or Schur vectors.
493 c See Remarks below.
494 c
495 c HOWMNY Character*1 (INPUT)
496 c Specifies the form of the basis for the invariant subspace
497 c corresponding to the converged Ritz values that is to be computed.
498 c
499 c = 'A': Compute NEV Ritz vectors;
500 c = 'P': Compute NEV Schur vectors;
501 c = 'S': compute some of the Ritz vectors, specified
502 c by the logical array SELECT.
503 c
504 c SELECT Logical array of dimension NCV. (INPUT)
505 c If HOWMNY = 'S', SELECT specifies the Ritz vectors to be
506 c computed. To select the Ritz vector corresponding to a
507 c Ritz value (DR(j), DI(j)), SELECT(j) must be set to .TRUE..
508 c If HOWMNY = 'A' or 'P', SELECT is used as internal workspace.
509 c
510 c DR Double precision array of dimension NEV+1. (OUTPUT)
511 c If IPARAM(7) = 1,2 or 3 and SIGMAI=0.0 then on exit: DR contains
512 c the real part of the Ritz approximations to the eigenvalues of
513 c A*z = lambda*B*z.
514 c If IPARAM(7) = 3, 4 and SIGMAI is not equal to zero, then on exit:
515 c DR contains the real part of the Ritz values of OP computed by
516 c DNAUPD . A further computation must be performed by the user
517 c to transform the Ritz values computed for OP by DNAUPD to those
518 c of the original system A*z = lambda*B*z. See remark 3 below.
519 c
520 c DI Double precision array of dimension NEV+1. (OUTPUT)
521 c On exit, DI contains the imaginary part of the Ritz value
522 c approximations to the eigenvalues of A*z = lambda*B*z associated
523 c with DR.
524 c
525 c NOTE: When Ritz values are complex, they will come in complex
526 c conjugate pairs. If eigenvectors are requested, the
527 c corresponding Ritz vectors will also come in conjugate
528 c pairs and the real and imaginary parts of these are
529 c represented in two consecutive columns of the array Z
530 c (see below).
531 c
532 c Z Double precision N by NEV+1 array if RVEC = .TRUE. and HOWMNY = 'A'. (OUTPUT)
533 c On exit, if RVEC = .TRUE. and HOWMNY = 'A', then the columns of
534 c Z represent approximate eigenvectors (Ritz vectors) corresponding
535 c to the NCONV=IPARAM(5) Ritz values for eigensystem
536 c A*z = lambda*B*z.
537 c
538 c The complex Ritz vector associated with the Ritz value
539 c with positive imaginary part is stored in two consecutive
540 c columns. The first column holds the real part of the Ritz
541 c vector and the second column holds the imaginary part. The
542 c Ritz vector associated with the Ritz value with negative
543 c imaginary part is simply the complex conjugate of the Ritz vector
544 c associated with the positive imaginary part.
545 c
546 c If RVEC = .FALSE. or HOWMNY = 'P', then Z is not referenced.
547 c
548 c NOTE: If if RVEC = .TRUE. and a Schur basis is not required,
549 c the array Z may be set equal to first NEV+1 columns of the Arnoldi
550 c basis array V computed by DNAUPD . In this case the Arnoldi basis
551 c will be destroyed and overwritten with the eigenvector basis.
552 c
553 c LDZ Integer. (INPUT)
554 c The leading dimension of the array Z. If Ritz vectors are
555 c desired, then LDZ >= max( 1, N ). In any case, LDZ >= 1.
556 c
557 c SIGMAR Double precision (INPUT)
558 c If IPARAM(7) = 3 or 4, represents the real part of the shift.
559 c Not referenced if IPARAM(7) = 1 or 2.
560 c
561 c SIGMAI Double precision (INPUT)
562 c If IPARAM(7) = 3 or 4, represents the imaginary part of the shift.
563 c Not referenced if IPARAM(7) = 1 or 2. See remark 3 below.
564 c
565 c WORKEV Double precision work array of dimension 3*NCV. (WORKSPACE)
566 c
567 c **** The remaining arguments MUST be the same as for the ****
568 c **** call to DNAUPD that was just completed. ****
569 c
570 c NOTE: The remaining arguments
571 c
572 c BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM, IPNTR,
573 c WORKD, WORKL, LWORKL, INFO
574 c
575 c must be passed directly to DNEUPD following the last call
576 c to DNAUPD . These arguments MUST NOT BE MODIFIED between
577 c the the last call to DNAUPD and the call to DNEUPD .
578 c
579 c Three of these parameters (V, WORKL, INFO) are also output parameters:
580 c
581 c V Double precision N by NCV array. (INPUT/OUTPUT)
582 c
583 c Upon INPUT: the NCV columns of V contain the Arnoldi basis
584 c vectors for OP as constructed by DNAUPD .
585 c
586 c Upon OUTPUT: If RVEC = .TRUE. the first NCONV=IPARAM(5) columns
587 c contain approximate Schur vectors that span the
588 c desired invariant subspace. See Remark 2 below.
589 c
590 c NOTE: If the array Z has been set equal to first NEV+1 columns
591 c of the array V and RVEC=.TRUE. and HOWMNY= 'A', then the
592 c Arnoldi basis held by V has been overwritten by the desired
593 c Ritz vectors. If a separate array Z has been passed then
594 c the first NCONV=IPARAM(5) columns of V will contain approximate
595 c Schur vectors that span the desired invariant subspace.
596 c
597 c WORKL Double precision work array of length LWORKL. (OUTPUT/WORKSPACE)
598 c WORKL(1:ncv*ncv+3*ncv) contains information obtained in
599 c dnaupd . They are not changed by dneupd .
600 c WORKL(ncv*ncv+3*ncv+1:3*ncv*ncv+6*ncv) holds the
601 c real and imaginary part of the untransformed Ritz values,
602 c the upper quasi-triangular matrix for H, and the
603 c associated matrix representation of the invariant subspace for H.
604 c
605 c Note: IPNTR(9:13) contains the pointer into WORKL for addresses
606 c of the above information computed by dneupd .
607 c -------------------------------------------------------------
608 c IPNTR(9): pointer to the real part of the NCV RITZ values of the
609 c original system.
610 c IPNTR(10): pointer to the imaginary part of the NCV RITZ values of
611 c the original system.
612 c IPNTR(11): pointer to the NCV corresponding error bounds.
613 c IPNTR(12): pointer to the NCV by NCV upper quasi-triangular
614 c Schur matrix for H.
615 c IPNTR(13): pointer to the NCV by NCV matrix of eigenvectors
616 c of the upper Hessenberg matrix H. Only referenced by
617 c dneupd if RVEC = .TRUE. See Remark 2 below.
618 c -------------------------------------------------------------
619 c
620 c INFO Integer. (OUTPUT)
621 c Error flag on output.
622 c
623 c = 0: Normal exit.
624 c
625 c = 1: The Schur form computed by LAPACK routine dlahqr
626 c could not be reordered by LAPACK routine dtrsen .
627 c Re-enter subroutine dneupd with IPARAM(5)=NCV and
628 c increase the size of the arrays DR and DI to have
629 c dimension at least dimension NCV and allocate at least NCV
630 c columns for Z. NOTE: Not necessary if Z and V share
631 c the same space. Please notify the authors if this error
632 c occurs.
633 c
634 c = -1: N must be positive.
635 c = -2: NEV must be positive.
636 c = -3: NCV-NEV >= 2 and less than or equal to N.
637 c = -5: WHICH must be one of 'LM', 'SM', 'LR', 'SR', 'LI', 'SI'
638 c = -6: BMAT must be one of 'I' or 'G'.
639 c = -7: Length of private work WORKL array is not sufficient.
640 c = -8: Error return from calculation of a real Schur form.
641 c Informational error from LAPACK routine dlahqr .
642 c = -9: Error return from calculation of eigenvectors.
643 c Informational error from LAPACK routine dtrevc .
644 c = -10: IPARAM(7) must be 1,2,3,4.
645 c = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatible.
646 c = -12: HOWMNY = 'S' not yet implemented
647 c = -13: HOWMNY must be one of 'A' or 'P' if RVEC = .true.
648 c = -14: DNAUPD did not find any eigenvalues to sufficient
649 c accuracy.
650 c = -15: DNEUPD got a different count of the number of converged
651 c Ritz values than DNAUPD got. This indicates the user
652 c probably made an error in passing data from DNAUPD to
653 c DNEUPD or that the data was modified before entering
654 c DNEUPD
655 c
656 c\BeginLib
657 c
658 c\References:
659 c 1. D.C. Sorensen, "Implicit Application of Polynomial Filters in
660 c a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992),
661 c pp 357-385.
662 c 2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly
663 c Restarted Arnoldi Iteration", Rice University Technical Report
664 c TR95-13, Department of Computational and Applied Mathematics.
665 c 3. B.N. Parlett & Y. Saad, "Complex Shift and Invert Strategies for
666 c Real Matrices", Linear Algebra and its Applications, vol 88/89,
667 c pp 575-595, (1987).
668 c
669 c\Routines called:
670 c ivout ARPACK utility routine that prints integers.
671 c dmout ARPACK utility routine that prints matrices
672 c dvout ARPACK utility routine that prints vectors.
673 c dgeqr2 LAPACK routine that computes the QR factorization of
674 c a matrix.
675 c dlacpy LAPACK matrix copy routine.
676 c dlahqr LAPACK routine to compute the real Schur form of an
677 c upper Hessenberg matrix.
678 c dlamch LAPACK routine that determines machine constants.
679 c dlapy2 LAPACK routine to compute sqrt(x**2+y**2) carefully.
680 c dlaset LAPACK matrix initialization routine.
681 c dorm2r LAPACK routine that applies an orthogonal matrix in
682 c factored form.
683 c dtrevc LAPACK routine to compute the eigenvectors of a matrix
684 c in upper quasi-triangular form.
685 c dtrsen LAPACK routine that re-orders the Schur form.
686 c dtrmm Level 3 BLAS matrix times an upper triangular matrix.
687 c dger Level 2 BLAS rank one update to a matrix.
688 c dcopy Level 1 BLAS that copies one vector to another .
689 c ddot Level 1 BLAS that computes the scalar product of two vectors.
690 c dnrm2 Level 1 BLAS that computes the norm of a vector.
691 c dscal Level 1 BLAS that scales a vector.
692 c
693 c\Remarks
694 c
695 c 1. Currently only HOWMNY = 'A' and 'P' are implemented.
696 c
697 c Let trans(X) denote the transpose of X.
698 c
699 c 2. Schur vectors are an orthogonal representation for the basis of
700 c Ritz vectors. Thus, their numerical properties are often superior.
701 c If RVEC = .TRUE. then the relationship
702 c A * V(:,1:IPARAM(5)) = V(:,1:IPARAM(5)) * T, and
703 c trans(V(:,1:IPARAM(5))) * V(:,1:IPARAM(5)) = I are approximately
704 c satisfied. Here T is the leading submatrix of order IPARAM(5) of the
705 c real upper quasi-triangular matrix stored workl(ipntr(12)). That is,
706 c T is block upper triangular with 1-by-1 and 2-by-2 diagonal blocks;
707 c each 2-by-2 diagonal block has its diagonal elements equal and its
708 c off-diagonal elements of opposite sign. Corresponding to each 2-by-2
709 c diagonal block is a complex conjugate pair of Ritz values. The real
710 c Ritz values are stored on the diagonal of T.
711 c
712 c 3. If IPARAM(7) = 3 or 4 and SIGMAI is not equal zero, then the user must
713 c form the IPARAM(5) Rayleigh quotients in order to transform the Ritz
714 c values computed by DNAUPD for OP to those of A*z = lambda*B*z.
715 c Set RVEC = .true. and HOWMNY = 'A', and
716 c compute
717 c trans(Z(:,I)) * A * Z(:,I) if DI(I) = 0.
718 c If DI(I) is not equal to zero and DI(I+1) = - D(I),
719 c then the desired real and imaginary parts of the Ritz value are
720 c trans(Z(:,I)) * A * Z(:,I) + trans(Z(:,I+1)) * A * Z(:,I+1),
721 c trans(Z(:,I)) * A * Z(:,I+1) - trans(Z(:,I+1)) * A * Z(:,I),
722 c respectively.
723 c Another possibility is to set RVEC = .true. and HOWMNY = 'P' and
724 c compute trans(V(:,1:IPARAM(5))) * A * V(:,1:IPARAM(5)) and then an upper
725 c quasi-triangular matrix of order IPARAM(5) is computed. See remark
726 c 2 above.
727 c
728 c\Authors
729 c Danny Sorensen Phuong Vu
730 c Richard Lehoucq CRPC / Rice University
731 c Chao Yang Houston, Texas
732 c Dept. of Computational &
733 c Applied Mathematics
734 c Rice University
735 c Houston, Texas
736 c
737 c\SCCS Information: @(#)
738 c FILE: neupd.F SID: 2.7 DATE OF SID: 09/20/00 RELEASE: 2
739 */
740
741 int
main(void)742 main(void)
743 {
744 // matrix
745 NaiveMatrixHandler mh(5);
746 mh(1, 1) = 2.;
747 mh(1, 2) = 1.;
748 mh(2, 1) = 1.;
749 mh(2, 2) = 2.;
750 mh(2, 3) = 1.;
751 mh(3, 2) = 1.;
752 mh(3, 3) = 2.;
753 mh(3, 4) = 1.;
754 mh(4, 3) = 1.;
755 mh(4, 4) = 2.;
756 mh(4, 5) = 1.;
757 mh(5, 4) = 1.;
758 mh(5, 5) = 2.;
759
760 // shift
761 doublereal SIGMA = 1.;
762 doublereal SIGMAR = 0.;
763 doublereal SIGMAI = 0.;
764
765 // arpack-related vars
766 integer IDO; // 0 at first iteration; then set by dnaupd
767 const char *BMAT; // 'I' for standard problem
768 integer N; // size of problem
769 const char *WHICH; // "SM" to request smallest eigenvalues
770 integer NEV; // number of eigenvalues
771 doublereal TOL; // -1 to use machine precision
772 std::vector<doublereal> RESID; // residual vector (ignored if IDO==0)
773 integer NCV; // number of vectors in subspace
774 std::vector<doublereal> V; // Schur basis
775 integer LDV; // leading dimension of V (==N!)
776 integer IPARAM[11] = { 0 };
777 integer IPNTR[14] = { 0 };
778 std::vector<doublereal> WORKD;
779 std::vector<doublereal> WORKL;
780 integer LWORKL;
781 integer INFO;
782
783 IDO = 0;
784 BMAT = "I";
785 N = mh.iGetNumRows();
786 WHICH = "SM";
787 NEV = 2;
788 TOL = 0.;
789 RESID.resize(N, 0.);
790 NCV = 4;
791 V.resize(N*NCV, 0.);
792 LDV = N;
793 IPARAM[0] = 1;
794 IPARAM[2] = 300;
795 IPARAM[3] = 1;
796 IPARAM[6] = 1; // mode...
797 WORKD.resize(3*N, 0.);
798 LWORKL = 3*NCV*NCV + 6*NCV;
799 WORKL.resize(LWORKL, 0.);
800 INFO = 0;
801
802 int cnt = 0;
803 do {
804 std::cout << "before:" << std::endl
805 << "IDO=" << IDO << std::endl
806 << "BMAT=" << BMAT << std::endl
807 << "N=" << N << std::endl
808 << "WHICH=" << WHICH << std::endl
809 << "NEV=" << NEV << std::endl
810 << "TOL=" << TOL << std::endl
811 << "&RESID[0]=" << &RESID[0] << std::endl
812 << "NCV=" << NCV << std::endl
813 << "&V[0]=" << &V[0] << std::endl
814 << "LDV=" << LDV << std::endl
815 << "IPARAM=" << IPARAM[0]
816 << "," << IPARAM[1]
817 << "," << IPARAM[2]
818 << "," << IPARAM[3]
819 << "," << IPARAM[4]
820 << "," << IPARAM[5]
821 << "," << IPARAM[6]
822 << "," << IPARAM[7]
823 << "," << IPARAM[8]
824 << "," << IPARAM[9]
825 << "," << IPARAM[10]
826 << std::endl
827 << "IPNTR=" << IPNTR[0]
828 << "," << IPNTR[1]
829 << "," << IPNTR[2]
830 << "," << IPNTR[3]
831 << "," << IPNTR[4]
832 << "," << IPNTR[5]
833 << "," << IPNTR[6]
834 << "," << IPNTR[7]
835 << "," << IPNTR[8]
836 << "," << IPNTR[9]
837 << "," << IPNTR[10]
838 << "," << IPNTR[11]
839 << "," << IPNTR[12]
840 << "," << IPNTR[13]
841 << std::endl
842 << "&WORKD[0]=" << &WORKD[0] << std::endl
843 << "&WORKL[0]=" << &WORKL[0] << std::endl
844 << "LWORKL=" << LWORKL << std::endl
845 << "INFO=" << INFO << std::endl
846 << std::endl;
847
848 __FC_DECL__(dnaupd)(&IDO, &BMAT[0], &N, &WHICH[0], &NEV,
849 &TOL, &RESID[0], &NCV, &V[0], &LDV, &IPARAM[0], &IPNTR[0],
850 &WORKD[0], &WORKL[0], &LWORKL, &INFO);
851
852 std::cout << "after:" << std::endl
853 << "IDO=" << IDO << std::endl
854 << "BMAT=" << BMAT << std::endl
855 << "N=" << N << std::endl
856 << "WHICH=" << WHICH << std::endl
857 << "NEV=" << NEV << std::endl
858 << "TOL=" << TOL << std::endl
859 << "&RESID[0]=" << &RESID[0] << std::endl
860 << "NCV=" << NCV << std::endl
861 << "&V[0]=" << &V[0] << std::endl
862 << "LDV=" << LDV << std::endl
863 << "IPARAM=" << IPARAM[0]
864 << "," << IPARAM[1]
865 << "," << IPARAM[2]
866 << "," << IPARAM[3]
867 << "," << IPARAM[4]
868 << "," << IPARAM[5]
869 << "," << IPARAM[6]
870 << "," << IPARAM[7]
871 << "," << IPARAM[8]
872 << "," << IPARAM[9]
873 << "," << IPARAM[10]
874 << std::endl
875 << "IPNTR=" << IPNTR[0]
876 << "," << IPNTR[1]
877 << "," << IPNTR[2]
878 << "," << IPNTR[3]
879 << "," << IPNTR[4]
880 << "," << IPNTR[5]
881 << "," << IPNTR[6]
882 << "," << IPNTR[7]
883 << "," << IPNTR[8]
884 << "," << IPNTR[9]
885 << "," << IPNTR[10]
886 << "," << IPNTR[11]
887 << "," << IPNTR[12]
888 << "," << IPNTR[13]
889 << std::endl
890 << "&WORKD[0]=" << &WORKD[0] << std::endl
891 << "&WORKL[0]=" << &WORKL[0] << std::endl
892 << "LWORKL=" << LWORKL << std::endl
893 << "INFO=" << INFO << std::endl
894 << std::endl;
895
896 std::cout << "cnt=" << cnt << ": IDO=" << IDO << ", INFO=" << INFO << std::endl;
897
898 // compute Y = OP*X
899 MyVectorHandler X(N, &WORKD[IPNTR[0] - 1]);
900 MyVectorHandler Y(N, &WORKD[IPNTR[1] - 1]);
901
902 Y = X;
903 Y *= -SIGMA;
904 mh.MatVecIncMul(Y, X);
905
906 std::cout << "X:" << std::endl << X << std::endl;
907 std::cout << "Y:" << std::endl << Y << std::endl;
908
909 cnt++;
910 } while (IDO == 1 || IDO == -1);
911
912 if (INFO < 0) {
913 std::cerr << "error" << std::endl;
914 return 1;
915 }
916
917 logical RVEC = true;
918 const char *HOWMNY = "A";
919 std::vector<logical> SELECT(NCV);
920 std::vector<doublereal> DR(NEV + 1);
921 std::vector<doublereal> DI(NEV + 1);
922 std::vector<doublereal> Z(N*(NEV + 1));
923 integer LDZ = N;
924 std::vector<doublereal> WORKEV(3*NCV);
925
926 __FC_DECL__(dneupd)(&RVEC, &HOWMNY[0], &SELECT[0], &DR[0], &DI[0],
927 &Z[0], &LDZ, &SIGMAR, &SIGMAI, &WORKEV[0],
928 &BMAT[0], &N, &WHICH[0], &NEV,
929 &TOL, &RESID[0], &NCV, &V[0], &LDV, &IPARAM[0], &IPNTR[0],
930 &WORKD[0], &WORKL[0], &LWORKL, &INFO);
931
932 bool bCmplx = false;
933 bool bFirst;
934 for (integer nv = 0; nv < NEV; nv++) {
935 std::cout << "value " << std::setw(8) << nv << ":"
936 << std::setw(16) << DR[nv] + SIGMA
937 << (DI[nv] >= 0. ? " + " : " - ")
938 << std::setw(16) << std::abs(DI[nv]) << " i" << std::endl;
939
940 std::cout << "vector " << std::setw(8) << nv << ":"
941 << std::endl;
942 if (!bCmplx) {
943 if (DI[nv] != 0.) {
944 bCmplx = true;
945 bFirst = true;
946 }
947 }
948
949 if (!bCmplx) {
950 for (integer idx = 0; idx < N; idx++) {
951 std::cout
952 << std::setw(8) << idx << ":"
953 << std::setw(16) << Z[N*nv + idx]
954 << std::endl;
955 }
956
957 } else if (bFirst) {
958 bFirst = false;
959
960 for (integer idx = 0; idx < N; idx++) {
961 doublereal d = Z[N*(nv + 1) + idx];
962 std::cout
963 << std::setw(8) << idx << ":"
964 << std::setw(16) << Z[N*nv + idx]
965 << (d > 0. ? " + " : " - " )
966 << std::setw(16) << std::abs(d) << " i"
967 << std::endl;
968 }
969
970 } else {
971 bCmplx = false;
972
973 for (integer idx = 0; idx < N; idx++) {
974 doublereal d = Z[N*nv + idx];
975 std::cout
976 << std::setw(8) << idx << ":"
977 << std::setw(16) << Z[N*(nv - 1) + idx]
978 << (d > 0. ? " + " : " - " )
979 << std::setw(16) << std::abs(d) << " i"
980 << std::endl;
981 }
982 }
983 }
984
985 return 0;
986 }
987