Home

Context Navigation

dsdrv6.f

Last change on this file was 2696, checked in by kanani, 6 years ago
Merge of branch palm4u into trunk
File size: 17.5 KB

Line
1	program dsdrv6
2	c
3	c Program to illustrate the idea of reverse communication
4	c in Cayley mode for a generalized symmetric eigenvalue
5	c problem. The following program uses the two LAPACK subroutines
6	c dgttrf.f and dgttrs.f to factor and solve a tridiagonal system of
7	c equations.
8	c
9	c We implement example six of ex-sym.doc in DOCUMENTS directory
10	c
11	c\Example-6
12	c ... Suppose we want to solve Ax = lambdaM*x in inverse mode,
13	c where A and M are obtained by the finite element of the
14	c 1-dimensional discrete Laplacian
15	c d^2u / dx^2
16	c on the interval [0,1] with zero Dirichlet boundary condition
17	c using piecewise linear elements.
18	c
19	c ... OP = (inv[A-sigmaM])(A+sigma*M) and B = M.
20	c
21	c ... Use mode 5 of DSAUPD.
22	c
23	c\BeginLib
24	c
25	c\References:
26	c 1. R.G. Grimes, J.G. Lewis and H.D. Simon, "A Shifted Block Lanczos
27	c Algorithm for Solving Sparse Symmetric Generalized Eigenproblems",
28	c SIAM J. Matr. Anal. Apps., January (1993).
29	c
30	c\Routines called:
31	c dsaupd ARPACK reverse communication interface routine.
32	c dseupd ARPACK routine that returns Ritz values and (optionally)
33	c Ritz vectors.
34	c dgttrf LAPACK tridiagonal factorization routine.
35	c dgttrs LAPACK tridiagonal solve routine.
36	c daxpy Level 1 BLAS that computes y <- alpha*x+y.
37	c dcopy Level 1 BLAS that copies one vector to another.
38	c dscal Level 1 BLAS that scales a vector by a scalar.
39	c dnrm2 Level 1 BLAS that computes the norm of a vector.
40	c av Matrix vector multiplication routine that computes A*x.
41	c mv Matrix vector multiplication routine that computes M*x.
42	c
43	c\Author
44	c Danny Sorensen
45	c Richard Lehoucq
46	c Chao Yang
47	c Dept. of Computational &
48	c Applied Mathematics
49	c Rice University
50	c Houston, Texas
51	c
52	c\SCCS Information: @(#)
53	c FILE: sdrv6.F SID: 2.5 DATE OF SID: 10/17/00 RELEASE: 2
54	c
55	c\Remarks
56	c 1. None
57	c
58	c\EndLib
59	c------------------------------------------------------------------------
60	c
61	c %-----------------------------%
62	c \| Define leading dimensions \|
63	c \| for all arrays. \|
64	c \| MAXN: Maximum dimension \|
65	c \| of the A allowed. \|
66	c \| MAXNEV: Maximum NEV allowed \|
67	c \| MAXNCV: Maximum NCV allowed \|
68	c %-----------------------------%
69	c
70	integer maxn, maxnev, maxncv, ldv
71	parameter (maxn=256, maxnev=10, maxncv=25,
72	& ldv=maxn)
73	c
74	c %--------------%
75	c \| Local Arrays \|
76	c %--------------%
77	c
78	Double precision
79	& v(ldv,maxncv), workl(maxncv*(maxncv+8)),
80	& workd(3*maxn), d(maxncv,2), resid(maxn),
81	& ad(maxn), adl(maxn), adu(maxn), adu2(maxn),
82	& temp(maxn), ax(maxn), mx(maxn)
83	logical select(maxncv)
84	integer iparam(11), ipntr(11), ipiv(maxn)
85	c
86	c %---------------%
87	c \| Local Scalars \|
88	c %---------------%
89	c
90	character bmat1, which2
91	integer ido, n, nev, ncv, lworkl, info, j, ierr,
92	& nconv, maxitr, ishfts, mode
93	logical rvec
94	Double precision
95	& sigma, r1, r2, tol, h
96	c
97	c %------------%
98	c \| Parameters \|
99	c %------------%
100	c
101	Double precision
102	& zero, one, two, four, six
103	parameter (zero = 0.0D+0, one = 1.0D+0,
104	& four = 4.0D+0, six = 6.0D+0,
105	& two = 2.0D+0 )
106	c
107	c %-----------------------------%
108	c \| BLAS & LAPACK routines used \|
109	c %-----------------------------%
110	c
111	Double precision
112	& dnrm2
113	external daxpy, dcopy, dscal, dnrm2, dgttrf, dgttrs
114	c
115	c %--------------------%
116	c \| Intrinsic function \|
117	c %--------------------%
118	c
119	intrinsic abs
120	c
121	c %-----------------------%
122	c \| Executable statements \|
123	c %-----------------------%
124	c
125	c %--------------------------------------------------%
126	c \| The number N is the dimension of the matrix. A \|
127	c \| generalized eigenvalue problem is solved (BMAT = \|
128	c \| 'G'.) NEV is the number of eigenvalues to be \|
129	c \| approximated. Since the Cayley mode is used, \|
130	c \| WHICH is set to 'LM'. The user can modify NEV, \|
131	c \| NCV, SIGMA to solve problems of different sizes, \|
132	c \| and to get different parts of the spectrum. \|
133	c \| However, The following conditions must be \|
134	c \| satisfied: \|
135	c \| N <= MAXN, \|
136	c \| NEV <= MAXNEV, \|
137	c \| NEV + 1 <= NCV <= MAXNCV \|
138	c %--------------------------------------------------%
139	c
140	n = 100
141	nev = 4
142	ncv = 20
143	if ( n .gt. maxn ) then
144	print *, ' ERROR with _SDRV6: N is greater than MAXN '
145	go to 9000
146	else if ( nev .gt. maxnev ) then
147	print *, ' ERROR with _SDRV6: NEV is greater than MAXNEV '
148	go to 9000
149	else if ( ncv .gt. maxncv ) then
150	print *, ' ERROR with _SDRV6: NCV is greater than MAXNCV '
151	go to 9000
152	end if
153	bmat = 'G'
154	which = 'LM'
155	sigma = 1.5D+2
156	c
157	c %--------------------------------------------------%
158	c \| The work array WORKL is used in DSAUPD as \|
159	c \| workspace. Its dimension LWORKL is set as \|
160	c \| illustrated below. The parameter TOL determines \|
161	c \| the stopping criterion. If TOL<=0, machine \|
162	c \| precision is used. The variable IDO is used for \|
163	c \| reverse communication and is initially set to 0. \|
164	c \| Setting INFO=0 indicates that a random vector is \|
165	c \| generated in DSAUPD to start the Arnoldi \|
166	c \| iteration. \|
167	c %--------------------------------------------------%
168	c
169	lworkl = ncv*(ncv+8)
170	tol = zero
171	ido = 0
172	info = 0
173	c
174	c %---------------------------------------------------%
175	c \| This program uses exact shifts with respect to \|
176	c \| the current Hessenberg matrix (IPARAM(1) = 1). \|
177	c \| IPARAM(3) specifies the maximum number of Arnoldi \|
178	c \| iterations allowed. Mode 5 specified in the \|
179	c \| documentation of DSAUPD is used (IPARAM(7) = 5). \|
180	c \| All these options may be changed by the user. For \|
181	c \| details, see the documentation in DSAUPD. \|
182	c %---------------------------------------------------%
183	c
184	ishfts = 1
185	maxitr = 300
186	mode = 5
187	c
188	iparam(1) = ishfts
189	iparam(3) = maxitr
190	iparam(7) = mode
191	c
192	c %------------------------------------------------------%
193	c \| Call LAPACK routine to factor (A-sigma*M). The \|
194	c \| stiffness matrix A is the 1-d discrete Laplacian. \|
195	c \| The mass matrix M is the associated mass matrix \|
196	c \| arising from using piecewise linear finite elements \|
197	c \| on the interval [0, 1]. \|
198	c %------------------------------------------------------%
199	c
200	h = one / dble(n+1)
201	r1 = (four / six) * h
202	r2 = (one / six) * h
203	do 20 j=1,n
204	ad(j) = two / h - sigma * r1
205	adl(j) = -one / h - sigma * r2
206	20 continue
207	call dcopy (n, adl, 1, adu, 1)
208	call dgttrf (n, adl, ad, adu, adu2, ipiv, ierr)
209	if (ierr .ne. 0) then
210	print *, ' '
211	print *, ' Error with _gttrf in _SDRV6.'
212	print *, ' '
213	go to 9000
214	end if
215	c
216	c %-------------------------------------------%
217	c \| M A I N L O O P (Reverse communication) \|
218	c %-------------------------------------------%
219	c
220	10 continue
221	c
222	c %---------------------------------------------%
223	c \| Repeatedly call the routine DSAUPD and take \|
224	c \| actions indicated by parameter IDO until \|
225	c \| either convergence is indicated or maxitr \|
226	c \| has been exceeded. \|
227	c %---------------------------------------------%
228	c
229	call dsaupd ( ido, bmat, n, which, nev, tol, resid, ncv,
230	& v, ldv, iparam, ipntr, workd, workl, lworkl,
231	& info )
232	c
233	if (ido .eq. -1) then
234	c
235	c %-------------------------------------------------------%
236	c \| Perform y <--- OPx = (inv[A-SIGMAM])(A+SIGMAM)*x \|
237	c \| to force starting vector into the range of OP. The \|
238	c \| user should provide his/her matrix vector (Ax, Mx) \|
239	c \| multiplication routines and a linear system solver \|
240	c \| here. The matrix vector multiplication routine takes \|
241	c \| workd(ipntr(1)) as the input vector. The final \|
242	c \| result is returned to workd(ipntr(2)). \|
243	c %-------------------------------------------------------%
244	c
245	call av (n, workd(ipntr(1)), workd(ipntr(2)))
246	call mv (n, workd(ipntr(1)), temp)
247	call daxpy(n, sigma, temp, 1, workd(ipntr(2)), 1)
248	c
249	call dgttrs ('Notranspose', n, 1, adl, ad, adu, adu2, ipiv,
250	& workd(ipntr(2)), n, ierr)
251	if (ierr .ne. 0) then
252	print *, ' '
253	print *, ' Error with _gttrs in _SDRV6.'
254	print *, ' '
255	go to 9000
256	end if
257	c
258	c %-----------------------------------------%
259	c \| L O O P B A C K to call DSAUPD again. \|
260	c %-----------------------------------------%
261	c
262	go to 10
263	c
264	else if (ido .eq. 1) then
265	c
266	c %----------------------------------------------------%
267	c \| Perform y <-- OPx = inv[A-SIGMAM](A+SIGMAM)*x. \|
268	c \| M*x has been saved in workd(ipntr(3)). The user \|
269	c \| need the matrix vector multiplication (A*x) \|
270	c \| routine and a linear system solver here. The \|
271	c \| matrix vector multiplication routine takes \|
272	c \| workd(ipntr(1)) as the input, and the result is \|
273	c \| combined with workd(ipntr(3)) to form the input \|
274	c \| for the linear system solver. The final result is \|
275	c \| returned to workd(ipntr(2)). \|
276	c %----------------------------------------------------%
277	c
278	call av (n, workd(ipntr(1)), workd(ipntr(2)))
279	call daxpy(n, sigma, workd(ipntr(3)), 1, workd(ipntr(2)), 1)
280	call dgttrs ('Notranspose', n, 1, adl, ad, adu, adu2, ipiv,
281	& workd(ipntr(2)), n, ierr)
282	if (ierr .ne. 0) then
283	print *, ' '
284	print *, ' Error with _gttrs in _SDRV6. '
285	print *, ' '
286	go to 9000
287	end if
288	c
289	c %-----------------------------------------%
290	c \| L O O P B A C K to call DSAUPD again. \|
291	c %-----------------------------------------%
292	c
293	go to 10
294	c
295	else if (ido .eq. 2) then
296	c
297	c %--------------------------------------------%
298	c \| Perform y <--- M*x. \|
299	c \| Need matrix vector multiplication routine \|
300	c \| here that takes workd(ipntr(1)) as input \|
301	c \| and returns the result to workd(ipntr(2)). \|
302	c %--------------------------------------------%
303	c
304	call mv (n, workd(ipntr(1)), workd(ipntr(2)))
305	c
306	c %-----------------------------------------%
307	c \| L O O P B A C K to call DSAUPD again. \|
308	c %-----------------------------------------%
309	c
310	go to 10
311	c
312	end if
313	c
314	c %-----------------------------------------%
315	c \| Either we have convergence, or there is \|
316	c \| an error. \|
317	c %-----------------------------------------%
318	c
319	if ( info .lt. 0 ) then
320	c
321	c %--------------------------%
322	c \| Error message, check the \|
323	c \| documentation in DSAUPD \|
324	c %--------------------------%
325	c
326	print *, ' '
327	print *, ' Error with _saupd, info = ',info
328	print *, ' Check the documentation of _saupd. '
329	print *, ' '
330	c
331	else
332	c
333	c %-------------------------------------------%
334	c \| No fatal errors occurred. \|
335	c \| Post-Process using DSEUPD. \|
336	c \| \|
337	c \| Computed eigenvalues may be extracted. \|
338	c \| \|
339	c \| Eigenvectors may also be computed now if \|
340	c \| desired. (indicated by rvec = .true.) \|
341	c %-------------------------------------------%
342	c
343	rvec = .true.
344	c
345	call dseupd ( rvec, 'All', select, d, v, ldv, sigma,
346	& bmat, n, which, nev, tol, resid, ncv, v, ldv,
347	& iparam, ipntr, workd, workl, lworkl, ierr )
348	c
349	c %----------------------------------------------%
350	c \| Eigenvalues are returned in the first column \|
351	c \| of the two dimensional array D and the \|
352	c \| corresponding eigenvectors are returned in \|
353	c \| the first NEV columns of the two dimensional \|
354	c \| array V if requested. Otherwise, an \|
355	c \| orthogonal basis for the invariant subspace \|
356	c \| corresponding to the eigenvalues in D is \|
357	c \| returned in V. \|
358	c %----------------------------------------------%
359	c
360	if ( ierr .ne. 0 ) then
361	c
362	c %------------------------------------%
363	c \| Error condition: \|
364	c \| Check the documentation of DSEUPD. \|
365	c %------------------------------------%
366	c
367	print *, ' '
368	print *, ' Error with _seupd, info = ', ierr
369	print *, ' Check the documentation of _seupd '
370	print *, ' '
371	c
372	else
373	c
374	c %---------------------------%
375	c \| Compute the residual norm \|
376	c \| \|
377	c \| \|\| Ax - lambdax \|\| \|
378	c \| \|
379	c \| for the NCONV accurately \|
380	c \| computed eigenvalues and \|
381	c \| eigenvectors. (iparam(5) \|
382	c \| indicates how many are \|
383	c \| accurate to the requested \|
384	c \| tolerance) \|
385	c %---------------------------%
386	c
387	nconv = iparam(5)
388	do 30 j=1, nconv
389	call av(n, v(1,j), ax)
390	call mv(n, v(1,j), mx)
391	call daxpy (n, -d(j,1), mx, 1, ax, 1)
392	d(j,2) = dnrm2(n, ax, 1)
393	d(j,2) = d(j,2) / abs(d(j,1))
394	30 continue
395	c
396	call dmout(6, nconv, 2, d, maxncv, -6,
397	& 'Ritz values and relative residuals')
398	c
399	end if
400	c
401	c %------------------------------------------%
402	c \| Print additional convergence information \|
403	c %------------------------------------------%
404	c
405	if ( info .eq. 1) then
406	print *, ' '
407	print *, ' Maximum number of iterations reached.'
408	print *, ' '
409	else if ( info .eq. 3) then
410	print *, ' '
411	print *, ' No shifts could be applied during implicit',
412	& ' Arnoldi update, try increasing NCV.'
413	print *, ' '
414	end if
415	c
416	print *, ' '
417	print *, ' _SDRV6 '
418	print *, ' ====== '
419	print *, ' '
420	print *, ' Size of the matrix is', n
421	print *, ' The number of Ritz values requested is', nev
422	print *, ' The number of Arnoldi vectors generated',
423	& ' (NCV) is ', ncv
424	print *, ' What portion of the spectrum: ', which
425	print *, ' The number of converged Ritz values is ',
426	& nconv
427	print *, ' The number of Implicit Arnoldi update',
428	& ' iterations taken is', iparam(3)
429	print , ' The number of OPx is ', iparam(9)
430	print *, ' The convergence criterion is ', tol
431	print *, ' '
432	c
433	end if
434	c
435	c %---------------------------%
436	c \| Done with program dsdrv6. \|
437	c %---------------------------%
438	c
439	9000 continue
440	c
441	end
442	c
443	c------------------------------------------------------------------------
444	c Matrix vector subroutine
445	c where the matrix used is the 1 dimensional mass matrix
446	c arising from using the piecewise linear finite element
447	c on the interval [0,1].
448	c
449	subroutine mv (n, v, w)
450	integer n, j
451	Double precision
452	& v(n), w(n), one, four, six, h
453	parameter (one = 1.0D+0, four = 4.0D+0,
454	& six = 6.0D+0)
455	c
456	w(1) = four*v(1) + v(2)
457	do 100 j = 2,n-1
458	w(j) = v(j-1) + four*v(j) + v(j+1)
459	100 continue
460	j = n
461	w(j) = v(j-1) + four*v(j)
462	c
463	c Scale the vector w by h.
464	c
465	h = one / (six*dble(n+1))
466	call dscal(n, h, w, 1)
467	return
468	end
469	c
470	c------------------------------------------------------------------------
471	c Matrix vector subroutine
472	c where the matrix is the stiffness matrix obtained from the
473	c finite element discretization of the 1-dimensional discrete Laplacian
474	c on the interval [0,1] with zero Dirichlet boundary condition
475	c using piecewise linear elements.
476	c
477	subroutine av (n, v, w)
478	integer n, j
479	Double precision
480	& v(n), w(n), one, two, h
481	parameter (one = 1.0D+0, two = 2.0D+0)
482	c
483	w(1) = two*v(1) - v(2)
484	do 100 j = 2,n-1
485	w(j) = - v(j-1) + two*v(j) - v(j+1)
486	100 continue
487	j = n
488	w(j) = - v(j-1) + two*v(j)
489	c
490	c Scale the vector w by (1/h).
491	c
492	h = one / dble(n+1)
493	call dscal(n, one/h, w, 1)
494	return
495	end

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/examples/dsdrv6.f

Download in other formats: