Context Navigation

pcg.F90 @ e083829

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

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

Download in other formats:

Original Format