Context Navigation

pcg.F90 @ ea465d3

Revision ea465d3, 21.7 KB checked in by Oleg Batrashev <ogbash@…>, 2 years ago (diff)
Move cdat and cdat setup under Coarse Preconditioner.
Property mode set to `100644`

Line
1	! DOUG - Domain decomposition On Unstructured Grids
2	! Copyright (C) 1998-2006 Faculty of Computer Science, University of Tartu and
3	! Department of Mathematics, University of Bath
4	!
5	! This library is free software; you can redistribute it and/or
6	! modify it under the terms of the GNU Lesser General Public
7	! License as published by the Free Software Foundation; either
8	! version 2.1 of the License, or (at your option) any later version.
9	!
10	! This library 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 GNU
13	! Lesser General Public License for more details.
14	!
15	! You should have received a copy of the GNU Lesser General Public
16	! License along with this library; if not, write to the Free Software
17	! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18	! or contact the authors (University of Tartu, Faculty of Computer Science, Chair
19	! of Distributed Systems, Liivi 2, 50409 Tartu, Estonia, http://dougdevel.org,
20	! mailto:info(at)dougdevel.org)
21
22	!!----------------------------------
23	!! Preconditioned Conjugate Gradient
24	!!----------------------------------
25	module pcg_mod
26
27	use ConvInf_mod
28	use SpMtx_mods
29	use Mesh_class
30	use globals
31	use subsolvers
32	use Preconditioner_mod
33
34	implicit none
35
36	#include<doug_config.h>
37
38	#ifdef D_COMPLEX
39	#define float complex
40	#else
41	#define float real
42	#endif
43
44	private :: &
45	preconditioner, &
46	msolve
47
48	! profiling
49	real(kind=rk), private, save :: time_preconditioner = 0
50
51	contains
52
53	subroutine prec2Level(prepare,A,CP,sol,rhs,res,CoarseMtx_,Restrict,isFirstIter)
54	use CoarseAllgathers
55
56	logical, intent(in) :: prepare
57	type(SpMtx) :: A !< sparse system matrix
58	type(CoarsePreconditioner),intent(inout) :: CP
59	real(kind=rk),dimension(:),pointer :: sol !< solution
60	real(kind=rk),dimension(:),pointer :: rhs !< right hand side
61	real(kind=rk),dimension(:),pointer :: res !< residual vector
62	type(SpMtx) :: CoarseMtx_ !< Coarse matrix
63	type(SpMtx) :: Restrict !< Restriction matrix
64	logical :: isFirstIter
65
66	real(kind=rk),dimension(:),pointer,save :: csol,crhs,clrhs,tmpsol2,tmpsol
67	integer :: ol
68
69	ol=max(sctls%overlap,sctls%smoothers)
70
71	if (prepare) then
72	if (sctls%verbose>4) write(stream,*) "Preparing 2. level"
73	call prec2Level_prepare()
74	else
75	if (sctls%verbose>4) write(stream,*) "Solving 2. level"
76	call prec2Level_solve()
77	end if
78
79	contains
80	subroutine prec2Level_exchangeMatrix()
81	integer :: ol
82	logical :: add
83
84	if (sctls%verbose>4) write(stream,*) "Exchanging coarse matrix"
85
86	ol=max(sctls%overlap,sctls%smoothers)
87	if (ol==0) then
88	add=.false.
89	else
90	add=.true.
91	endif
92
93	call AllSendCoarseMtx(CP%cdat%LAC,CoarseMtx_,CP%cdat%lg_cfmap,&
94	CP%cdat%ngfc,CP%cdat%nprocs,CP%cdat%send)
95	call AllRecvCoarseMtx(CoarseMtx_,CP%cdat%send,add=add) ! Recieve it
96
97	end subroutine prec2Level_exchangeMatrix
98
99	subroutine prec2Level_prepare()
100	if (isFirstIter) then
101	if (CP%cdat%active) then
102	! First iteration - send matrix
103	call prec2Level_exchangeMatrix()
104	end if
105
106	! after coarse matrix is received allocate coarse vectors
107	if (associated(crhs)) then
108	if (size(crhs)/=CoarseMtx_%ncols) then
109	deallocate(csol)
110	deallocate(crhs)
111	allocate(crhs(CoarseMtx_%ncols))
112	allocate(csol(CoarseMtx_%nrows))
113	endif
114	else
115	allocate(crhs(CoarseMtx_%ncols))
116	allocate(csol(CoarseMtx_%nrows))
117	endif
118	! allocate memory for vector
119	if (CP%cdat%active) then
120	allocate(clrhs(CP%cdat%nlfc))
121	else
122	allocate(clrhs(CP%cdat_vec%nlfc))
123	end if
124	end if
125
126	if (.not.associated(tmpsol)) then
127	!allocate(tmpsol(A%nrows))
128	allocate(tmpsol(size(rhs)))
129	endif
130
131	if (CP%cdat%active) then
132	if (sctls%verbose>6) write(stream,*) "Restricting into local coarse vector", size(clrhs)
133	! Send coarse vector
134	call SpMtx_Ax(clrhs,Restrict,rhs,dozero=.true.) ! restrict <RA>
135	if (CP%cdat_vec%active) then
136	call AllSendCoarseVector(clrhs,CP%cdat_vec%nprocs,CP%cdat_vec%cdisps,&
137	CP%cdat_vec%send,useprev=.not.isFirstIter)
138	else
139	call AllSendCoarseVector(clrhs,CP%cdat%nprocs,CP%cdat%cdisps,&
140	CP%cdat%send,useprev=.not.isFirstIter)
141	endif
142	end if ! CP%cdat%active
143	end subroutine prec2Level_prepare
144
145	subroutine prec2Level_solve()
146	if (.not.CP%cdat%active) then ! 1 processor case
147	if (sctls%method>1.and.sctls%method/=5) then ! multiplicative Schwarz
148	call SpMtx_Ax(crhs,Restrict,res,dozero=.true.) ! restriction
149	else
150	call SpMtx_Ax(crhs,Restrict,rhs,dozero=.true.) ! restriction
151	endif
152	end if
153
154	!csol=0.0_rk
155	if (CP%cdat%active) then
156	! Recieve the vector for solve
157	if (CP%cdat_vec%active) then
158	call AllRecvCoarseVector(crhs,CP%cdat_vec%nprocs,&
159	CP%cdat_vec%cdisps,CP%cdat_vec%glg_cfmap,CP%cdat_vec%send)
160	else
161	call AllRecvCoarseVector(crhs,CP%cdat%nprocs,&
162	CP%cdat%cdisps,CP%cdat%glg_cfmap,CP%cdat%send)
163	endif
164	!call MPI_BARRIER(MPI_COMM_WORLD,i)
165	end if
166
167	if (isFirstIter) then
168	write (stream,*) &
169	'factorising coarse matrix of size',CoarseMtx_%nrows, &
170	' and nnz:',CoarseMtx_%nnz
171	CoarseMtx_%subsolve_id=0
172	end if
173
174	! Coarse solve
175	call sparse_singlesolve(CoarseMtx_%subsolve_id,csol,crhs,&
176	nfreds=CoarseMtx_%nrows, &
177	nnz=CoarseMtx_%nnz, &
178	indi=CoarseMtx_%indi, &
179	indj=CoarseMtx_%indj, &
180	val=CoarseMtx_%val)
181	if (isFirstIter) then
182	CoarseMtx_%indi=CoarseMtx_%indi+1
183	CoarseMtx_%indj=CoarseMtx_%indj+1
184	end if
185
186	if (CP%cdat_vec%active) then
187	call Vect_remap(csol,clrhs,CP%cdat%gl_cfmap,dozero=.true.)
188	call SpMtx_Ax(tmpsol,Restrict,clrhs,dozero=.true.,transp=.true.)
189
190	elseif (CP%cdat%active) then
191	call Vect_remap(csol,clrhs,CP%cdat%gl_cfmap,dozero=.true.)
192	call SpMtx_Ax(tmpsol,Restrict,clrhs,dozero=.true.,transp=.true.)
193	else
194	call SpMtx_Ax(tmpsol,Restrict,csol,dozero=.true.,transp=.true.) ! interpolation
195	endif
196
197	if (sctls%method==1) then
198	!call Print_Glob_Vect(sol,M,'sol===',chk_endind=M%ninner)
199	sol=sol+tmpsol
200	elseif (sctls%method==3) then ! fully multiplicative Schwarz
201	sol(1:A%nrows)=sol(1:A%nrows)+tmpsol(1:A%nrows)
202	! calculate the residual:
203	call SpMtx_Ax(res,A,sol,dozero=.true.) !
204	res=rhs-res
205	elseif (sctls%method==2.and.ol>0) then ! fully multiplicative Schwarz
206	sol(1:A%nrows)=tmpsol(1:A%nrows)
207	! calculate the residual:
208	call SpMtx_Ax(res,A,sol,dozero=.true.) !
209	res=rhs-res
210	endif
211	if (((ol==0.and.sctls%method==2).or.sctls%method==5).and.sctls%levels>1) then
212	! multiplicative on fine level, additive with coarse level:
213	sol(1:A%nrows)=sol(1:A%nrows)+tmpsol(1:A%nrows)
214	endif
215	end subroutine prec2Level_solve
216
217	end subroutine prec2Level
218
219	!-------------------------------
220	!> Make preconditioner
221	!-------------------------------
222	subroutine preconditioner(sol,A,rhs,M,finePrec,coarsePrec,&
223	A_ghost,CoarseMtx_,Restrict,refactor,bugtrack_)
224	use CoarseAllgathers
225	use CoarseMtx_mod
226	use Vect_mod
227	implicit none
228	real(kind=rk),dimension(:),pointer :: sol !< solution
229	type(SpMtx) :: A !< sparse system matrix
230	real(kind=rk),dimension(:),pointer :: rhs !< right hand side
231	type(Mesh),intent(in) :: M !< Mesh
232	type(FinePreconditioner),intent(inout) :: finePrec !< fine level preconditioner
233	type(CoarsePreconditioner),intent(inout) :: coarsePrec !< coarse level preconditioner
234	real(kind=rk),dimension(:),pointer :: res !< residual vector, allocated
235	!! here for multiplicative Schwarz
236	type(SpMtx),optional :: A_ghost !< matr@interf.
237	type(SpMtx),optional :: CoarseMtx_ !< Coarse matrix
238	type(SpMtx),optional :: Restrict !< Restriction matrix
239	logical,intent(inout),optional :: refactor
240	logical,optional :: bugtrack_
241	! ----- local: ------
242	integer :: i
243	logical :: add,bugtrack
244	real(kind=rk) :: t1
245	logical :: isFirstIter
246
247	if (sctls%verbose>4) write(stream,*) "Applying preconditioner"
248	t1 = MPI_WTime()
249
250	if (present(bugtrack_)) then
251	bugtrack=bugtrack_
252	else
253	bugtrack=.false.
254	endif
255	! ----------------------------
256	isFirstIter = .false.
257	if (present(refactor)) isFirstIter = refactor
258
259	if (sctls%method==0) then
260	sol=rhs
261	return
262	endif
263
264	sol=0.0_rk
265	if (sctls%method>1) then ! For multiplicative Schwarz method...:
266	allocate(res(size(rhs)))
267	endif
268
269	if (sctls%levels>1) then
270	call prec2Level(.true.,A,coarsePrec,sol,rhs,res,CoarseMtx_,Restrict,isFirstIter)
271	end if
272
273	! first level prec
274	call FinePreconditioner_apply(finePrec,sol,rhs)
275
276	if (sctls%levels>1) then
277	call prec2Level(.false.,A,coarsePrec,sol,rhs,res,CoarseMtx_,Restrict,isFirstIter)
278	end if
279
280	time_preconditioner = time_preconditioner + MPI_WTime()-t1
281
282	end subroutine preconditioner
283
284	!--------------------------
285	! Solve
286	!--------------------------
287	subroutine msolve (m, r, z)
288	implicit none
289	real(kind=rk), dimension(:), intent(in) :: m, r
290	real(kind=rk), dimension(:), intent(in out) :: z
291	integer :: i
292
293	do i = 1, size(r)
294	z(i) = r(i) !/ m(i)
295	end do
296	end subroutine msolve
297
298	!--------------------------
299	!> Preconditioned conjugent gradient method with eigenvalues
300	!--------------------------
301	subroutine pcg_weigs (A,b,x,Msh,finePrec,coarsePrec,it,cond_num,A_interf_,tol_,maxit_, &
302	x0_,solinf,resvects_,CoarseMtx_,Restrict,refactor_)
303	use CoarseAllgathers
304
305	implicit none
306
307	type(SpMtx),intent(in out) :: A !< System matrix (sparse)
308	float(kind=rk),dimension(:),pointer :: b !< right hand side
309	float(kind=rk),dimension(:),pointer :: x !< Solution
310	type(Mesh),intent(in) :: Msh !< Mesh - aux data for Ax operation
311	type(FinePreconditioner),intent(inout) :: finePrec !< fine level preconditioner
312	type(CoarsePreconditioner),intent(inout) :: coarsePrec !< coarse level preconditioner
313
314	integer,intent(out) :: it
315	real(kind=rk),intent(out) :: cond_num
316
317	! Optional arguments
318	type(SpMtx),intent(in out),optional :: A_interf_ !< matr@interf.
319	real(kind=rk), intent(in), optional :: tol_ !< Tolerance of the method
320	integer, intent(in), optional :: maxit_ !< Max number of iterations
321	float(kind=rk), dimension(:), intent(in), optional :: x0_ !< Initial guess
322	type(ConvInf), intent(in out), optional :: solinf !< Solution statistics
323	logical, intent(in), optional :: resvects_ !< Fill in the 'resvect' or not
324	type(SpMtx),optional :: CoarseMtx_ !< Coarse matrix
325	type(SpMtx),optional :: Restrict !< Restriction mtx
326	logical,intent(in),optional :: refactor_
327
328	! Local variables
329	logical :: refactor
330
331	real(kind=rk) :: tol ! Tolerance
332	integer :: maxit ! Max number of iterations
333	logical :: resvects=.false.
334
335	real(kind=rk),dimension(:),pointer,save :: r, p, q, m, z, b_k
336	real(kind=rk),dimension(:),pointer,save :: ztmp !TODO remove me
337	real(kind=rk) :: rho_curr, rho_prev, init_norm, res_norm
338	real(kind=rk) :: tmp,ratio_norm
339	integer :: nfreds
340	integer :: ierr
341	real(kind=rk),dimension(:),pointer,save :: dd,ee,alpha,beta
342	!logical,parameter :: bugtrack=.true.
343	logical,parameter :: bugtrack=.false.
344	! profiling
345	real(8) :: time_iterations, t1
346
347	write(stream,'(/a)') 'Preconditioned conjugate gradient:'
348	if (present(refactor_).and.refactor_) then
349	refactor=.true.
350	else
351	refactor=.false.
352	endif
353
354	if (size(b) /= size(x)) &
355	call DOUG_abort('[pcg] : SEVERE : size(b) /= size(x)',-1)
356
357	nfreds=size(b)
358	if (.not.associated(r)) then
359	allocate(r(nfreds))
360	allocate(p(nfreds))
361	allocate(q(nfreds))
362	allocate(m(nfreds))
363	allocate(z(nfreds))
364	allocate(b_k(nfreds))
365	endif
366
367	tol = sctls%solve_tolerance
368	if (present(tol_)) &
369	tol = tol_
370
371	maxit = sctls%solve_maxiters
372	if (maxit==-1) maxit=100000
373	if (present(maxit_)) &
374	maxit = maxit_
375
376	if (present(resvects_).and.(.not.present(solinf))) &
377	call DOUG_abort('[pcg] : SEVERE : "resvects_" must be given along'//&
378	' with "solinf".',-1)
379
380	! Allocate convergence info structure
381	if (present(solinf).and.(.not.present(resvects_))) then
382	call ConvInf_Init(solinf)
383	else if (present(solinf).and.&
384	(present(resvects_).and.(resvects_.eqv.(.true.)))) then
385	call ConvInf_Init(solinf, maxit)
386	resvects = .true.
387	end if
388
389
390	! Initialise auxiliary data structures
391	! to assist with pmvm
392	call pmvmCommStructs_init(A, Msh)
393
394
395	if (bugtrack)call Print_Glob_Vect(x,Msh,'global x===')
396	call SpMtx_pmvm(r,A,x,Msh)
397
398	r = b - r
399	init_norm = Vect_dot_product(r,r)
400	if (init_norm == 0.0) &
401	init_norm = 1.0
402
403
404	write(stream,'(a,i6,a,e8.2)') 'maxit = ',maxit,', tol = ',tol
405	write(stream,'(a,e22.15)') 'init_norm = ', dsqrt(init_norm)
406
407	it = 0
408	ratio_norm = 1.0_rk
409	! arrays for eigenvalue calculation:
410	if (.not.associated(dd)) then
411	allocate(dd(maxit))
412	allocate(ee(maxit))
413	allocate(alpha(maxit))
414	allocate(beta(maxit))
415	endif
416
417	time_iterations = 0.
418	! iterations
419	do while((ratio_norm > tol*tol).and.(it < maxit))
420	it = it + 1
421	t1 = MPI_WTime()
422
423	if (bugtrack)call Print_Glob_Vect(r,Msh,'global r===',chk_endind=Msh%ninner)
424	call preconditioner(sol=z, &
425	A=A, &
426	rhs=r, &
427	M=Msh, &
428	finePrec=finePrec, &
429	coarsePrec=coarsePrec, &
430	A_ghost=A_interf_, &
431	CoarseMtx_=CoarseMtx_, &
432	Restrict=Restrict, &
433	refactor=refactor, &
434	bugtrack_=bugtrack)
435
436	refactor=.false.
437	if (bugtrack)call Print_Glob_Vect(z,Msh,'global bef comm z===',chk_endind=Msh%ninner)
438	if (sctls%method/=0) then
439	call Add_common_interf(z,A,Msh)
440	endif
441
442	!call Print_Glob_Vect(z,Msh,'global aft comm z===',chk_endind=Msh%ninner)
443	!call Print_Glob_Vect(z,Msh,'global z===')
444	! compute current rho
445	rho_curr = Vect_dot_product(r,z)
446
447	if (it == 1) then
448	p = z
449	else
450	beta(it) = rho_curr / rho_prev
451	p = z + beta(it) * p
452	end if
453	call SpMtx_pmvm(q,A, p, Msh)
454	if (bugtrack)call Print_Glob_Vect(q,Msh,'global q===')
455	! compute alpha
456	alpha(it) = rho_curr / Vect_dot_product(p,q)
457	x = x + alpha(it) * p
458	if (bugtrack)call Print_Glob_Vect(x,Msh,'global x===')
459	r = r - alpha(it) * q
460	if (bugtrack)call Print_Glob_Vect(r,Msh,'global r===')
461	rho_prev = rho_curr
462	! check
463	res_norm = Vect_dot_product(r,r)
464	!write (stream,*) "Norm is", res_norm
465	ratio_norm = res_norm / init_norm
466
467	time_iterations = time_iterations + MPI_WTime()-t1
468
469	if (ismaster()) &
470	write(stream, '(i5,a,e22.15)') it,': res_norm=',dsqrt(res_norm)
471	end do
472
473	if(pstream/=0) then
474	write(pstream, "(I0,':pcg iterations:',I0)") myrank, it
475	write(pstream, "(I0,':iterations time:',F0.3)") myrank, time_iterations
476	write(pstream, "(I0,':preconditioner time:',F0.3)") myrank, time_preconditioner
477	end if
478
479	if (ismaster()) then
480	call CalculateEigenvalues(it,dd,ee,alpha,beta,maxit)
481	cond_num=dd(it)/dd(1)
482	write(stream,'(a,i3,a,e10.4)') '#it:',it,' Cond#: ',cond_num
483	endif
484	!deallocate(beta)
485	!deallocate(alpha)
486	!deallocate(ee)
487	!deallocate(dd)
488	! Deallocate auxiliary data structures
489	! helped to assist with pmvm
490	! call pmvmCommStructs_destroy()
491
492	!deallocate(b_k)
493	!deallocate(z)
494	!deallocate(m)
495	!deallocate(q)
496	!deallocate(p)
497	!deallocate(r)
498	end subroutine pcg_weigs
499
500
501	subroutine CalculateEigenvalues(it,dd,ee,alpha,beta,MaxIt)
502	! Finds eigenvalues using Lanczos connection
503	implicit none
504	integer :: it,MaxIt
505	real(kind=rk),dimension(:),pointer :: dd,ee,alpha,beta
506	integer :: i,j,ierr
507
508	dd(1) = 1.0_rk/alpha(1)
509	ee(1) = 0.0_rk
510	do i=2,it
511	if (alpha(i)==0.or.beta(i)<0) then
512	print *,'Unable to compute eigenvalue number',i,' and '
513	do j=1,i
514	print *,j,' alpha:',alpha(j),' beta:',beta
515	enddo
516	return
517	else
518	dd(i) = 1.0_rk/alpha(i) + beta(i)/alpha(i-1)
519	ee(i) = -dsqrt(beta(i))/alpha(i-1)
520	endif
521	enddo
522	!call tql1(it,dd,ee,MaxIt,ierr)
523	call tql1(it,dd,ee,ierr)
524	if (ierr > 0) then
525	print *,'Cancelled after EigMaxIt while calculating eig',ierr
526	endif
527	end subroutine CalculateEigenvalues
528
529	subroutine tql1(n,d,e,ierr)
530	implicit none
531	integer :: i,j,l,m,n,ii,l1,l2,mml,ierr
532	real(kind=rk) :: d(n),e(n)
533	real(kind=rk) :: c,c2,c3,dl1,el1,f,g,h,p,r,s,s2,tst1,tst2 !,pythag
534	!
535	! this subroutine is a translation of the algol procedure tql1,
536	! num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and
537	! wilkinson.
538	! handbook for auto. comp., vol.ii-linear algebra, 227-240(1971).
539	!
540	! this subroutine finds the eigenvalues of a symmetric
541	! tridiagonal matrix by the ql method.
542	!
543	! on input
544	! n is the order of the matrix.
545	! d contains the diagonal elements of the input matrix.
546	! e contains the subdiagonal elements of the input matrix
547	! in its last n-1 positions. e(1) is arbitrary.
548	!
549	! on output
550	! d contains the eigenvalues in ascending order. if an
551	! error exit is made, the eigenvalues are correct and
552	! ordered for indices 1,2,...ierr-1, but may not be
553	! the smallest eigenvalues.
554	! e has been destroyed.
555	! ierr is set to
556	! zero for normal return,
557	! j if the j-th eigenvalue has not been
558	! determined after 30 iterations.
559	! calls pythag for sqrt(aa + bb) .
560	!
561	! questions and comments should be directed to burton s. garbow,
562	! mathematics and computer science div, argonne national laboratory
563	!
564	! this version dated august 1983.
565	!
566	! ------------------------------------------------------------------
567	!
568	ierr = 0
569	if (n .eq. 1) go to 1001
570	!
571	do 100 i = 2, n
572	100 e(i-1) = e(i)
573	!
574	f = 0.0e0
575	tst1 = 0.0e0
576	e(n) = 0.0e0
577
578	do 290 l = 1, n
579	j = 0
580	h = abs(d(l)) + abs(e(l))
581	if (tst1 .lt. h) tst1 = h
582	! .......... look for small sub-diagonal element ..........
583	do 110 m = l, n
584	tst2 = tst1 + abs(e(m))
585	if (tst2 .eq. tst1) go to 120
586	! .......... e(n) is always zero, so there is no exit
587	! throu gh the bottom of the loop ..........
588	110 continue
589
590	120 if (m .eq. l) go to 210
591	130 if (j .eq. 30 ) go to 1000
592	j = j + 1
593	! .......... form shift ..........
594	l1 = l + 1
595	l2 = l1 + 1
596	g = d(l)
597	p = (d(l1) - g) / (2.0e0 * e(l))
598	r = pythag(p,1.0e0_rk)
599	d(l) = e(l) / (p + sign(r,p))
600	d(l1) = e(l) * (p + sign(r,p))
601	dl1 = d(l1)
602	h = g - d(l)
603	if (l2 .gt. n) go to 145
604
605	do 140 i = l2, n
606	140 d(i) = d(i) - h
607
608	145 f = f + h
609	! .......... ql transformation ..........
610	p = d(m)
611	c = 1.0e0
612	c2 = c
613	el1 = e(l1)
614	s = 0.0e0
615	mml = m - l
616	! .......... for i=m-1 step -1 until l do -- ..........
617	do 200 ii = 1, mml
618	c3 = c2
619	c2 = c
620	s2 = s
621	i = m - ii
622	g = c * e(i)
623	h = c * p
624	r = pythag(p,e(i))
625	e(i+1) = s * r
626	s = e(i) / r
627	c = p / r
628	p = c * d(i) - s * g
629	d(i+1) = h + s * (c * g + s * d(i))
630	200 continue
631
632	p = -s * s2 * c3 * el1 * e(l) / dl1
633	e(l) = s * p
634	d(l) = c * p
635	tst2 = tst1 + abs(e(l))
636	if (tst2 .gt. tst1) go to 130
637	210 p = d(l) + f
638	! .......... order eigenvalues ..........
639	if (l .eq. 1) go to 250
640	! .......... for i=l step -1 until 2 do -- ..........
641	do 230 ii = 2, l
642	i = l + 2 - ii
643	if (p .ge. d(i-1)) go to 270
644	d(i) = d(i-1)
645	230 continue
646
647	250 i = 1
648	270 d(i) = p
649	290 continue
650
651	go to 1001
652	! .......... set error -- no convergence to an
653	! eigenvalue after 30 iterations ..........
654	1000 ierr = l
655	1001 return
656	end subroutine tql1
657
658	function pythag(a,b) result(pyt)
659	! finds dsqrt(a2+b2) without overflow or destructive underflow
660	implicit none
661	real(kind=rk) :: a,b,pyt
662	real(kind=rk) :: p,r,s,t,u
663
664	p = dmax1(dabs(a),dabs(b))
665	if (p .eq. 0.0d0) go to 20
666	r = (dmin1(dabs(a),dabs(b))/p)**2
667	10 continue
668	t = 4.0d0 + r
669	if (t .eq. 4.0d0) go to 20
670	s = r/t
671	u = 1.0d0 + 2.0d0*s
672	p = u*p
673	r = (s/u)*2 r
674	go to 10
675	20 pyt = p
676	return
677	end function pythag
678
679
680	end module pcg_mod
681
682

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

Download in other formats:

Original Format