Context Navigation

pcg.F90 @ 686ede6

Revision 686ede6, 29.7 KB checked in by Oleg Batrashev <ogbash@…>, 2 years ago (diff)
Clone commit
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 subsolvers_mult
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	integer,save :: coarseid=0
45	private :: &
46	preconditioner, &
47	msolve
48
49	! profiling
50	real(kind=rk), private, save :: time_preconditioner = 0
51
52	contains
53
54	!--------------------------------------------
55	! Preconditioned conjugate gradient method
56	!--------------------------------------------
57	subroutine pcg (A, b, x, Msh, tol_, maxit_, &
58	x0_, solinf, resvects_, CoarseMtx_)
59	implicit none
60
61	type(SpMtx),intent(in out) :: A ! System matrix (sparse)
62	float(kind=rk),dimension(:),pointer :: b ! RHS
63	float(kind=rk),dimension(:),pointer :: x ! Solution
64	! Mesh - aux data for Ax operation
65	type(Mesh), intent(in) :: Msh
66
67	! Optional arguments
68	! Tolerance of the method
69	real(kind=rk), intent(in), optional :: tol_
70	! Max number of iterations
71	integer, intent(in), optional :: maxit_
72	! Initial guess
73	float(kind=rk), dimension(:), intent(in), optional :: x0_
74	! Solution statistics
75	type(ConvInf), intent(in out), optional :: solinf
76	! Fill in the 'resvect' or not
77	logical, intent(in), optional :: resvects_
78	type(SpMtx),optional :: CoarseMtx_ ! Coarse matrix
79
80	real(kind=rk) :: tol ! Tolerance
81	integer :: maxit ! Max number of iterations
82	logical :: resvects=.false.
83
84	!real(kind=rk), dimension(size(b)) :: r, p, q, m, z, b_k
85	real(kind=rk),dimension(:),pointer :: r, p, q, m, z, b_k
86	real(kind=rk) :: rho_curr, rho_prev, init_norm, res_norm
87	real(kind=rk) :: alpha, beta, tmp, ratio_norm
88	integer :: iter,nfreds
89	integer :: ierr
90	! + test
91	float(kind=rk), dimension(:), pointer :: arr_copy, gv_tmp
92
93	!!$ real(kind=rk) :: r_max
94
95	write(stream,'(/a)') 'Preconditioned conjugate gradient:'
96
97	if (size(b) /= size(x)) &
98	call DOUG_abort('[pcg] : SEVERE : size(b) /= size(x)',-1)
99
100	!!$ if (present(x0_))
101	!!$ ! use method with initial guess
102
103
104	nfreds=size(b)
105	allocate(r(nfreds))
106	allocate(p(nfreds))
107	allocate(q(nfreds))
108	allocate(m(nfreds))
109	allocate(z(nfreds))
110	allocate(b_k(nfreds))
111
112
113	tol = sctls%solve_tolerance
114	if (present(tol_)) &
115	tol = tol_
116
117	maxit = sctls%solve_maxiters
118	if (maxit==-1) maxit=100000
119	if (present(maxit_)) &
120	maxit = maxit_
121
122	if (present(resvects_).and.(.not.present(solinf))) &
123	call DOUG_abort('[pcg] : SEVERE : "resvects_" must be given along'//&
124	' with "solinf".',-1)
125
126	! Allocate convergence info structure
127	if (present(solinf).and.(.not.present(resvects_))) then
128	call ConvInf_Init(solinf)
129	else if (present(solinf)) then
130	if (present(resvects_).and.(resvects_.eqv.(.true.))) then
131	call ConvInf_Init(solinf, maxit)
132	resvects = .true.
133	end if
134	end if
135
136	! Initialise auxiliary data structures
137	! to assist with pmvm
138	call pmvmCommStructs_init(A, Msh)
139
140	! + test
141	if (ismaster()) &
142	allocate(gv_tmp(Msh%ngf))
143	allocate(arr_copy(Msh%nlf))
144
145	! Build preconditioner
146	! call preconditioner(A, m, CoarseMtx_)
147	!!$ ! + test
148	!!$ arr_copy = m
149	!!$ call Vect_newToOldPerm(arr_copy)
150	!!$ gv_tmp = 0.0_rk
151	!!$ call Vect_Gather(arr_copy, gv_tmp, Msh)
152	!!$ if (ismaster()) &
153	!!$ call Vect_Print(gv_tmp, 'm ::')
154
155	!!$ ! + test
156	!!$ !b = 1.0_rk
157	!!$ arr_copy = b
158	!!$ call Vect_newToOldPerm(arr_copy)
159	!!$ call Vect_Print(arr_copy, 'b local (original ordering) ::')
160	!!$ gv_tmp = 0.0_rk
161	!!$ call Vect_Gather(arr_copy, gv_tmp, Msh)
162	!!$ if (ismaster()) &
163	!!$ call Vect_Print(gv_tmp, 'b global (original ordering) ::')
164
165	!!$ ! + test
166	!!$ arr_copy = x
167	!!$ call Vect_newToOldPerm(arr_copy)
168	!!$ call Vect_Print(arr_copy, 'x local (original ordering) ::')
169	!!$ gv_tmp = 0.0_rk
170	!!$ call Vect_Gather(arr_copy, gv_tmp, Msh)
171	!!$ if (ismaster()) &
172	!!$ call Vect_Print(gv_tmp, 'x global (original ordering) ::')
173
174	!r = b - SpMtx_pmvm(A, x, Msh)
175	call SpMtx_pmvm(r,A,x,Msh)
176	r = b - r
177	!call Vect_Print(r,'initial residual')
178
179	! + test
180	!!$ arr_copy = r
181	!!$ call Vect_newToOldPerm(arr_copy)
182	!!$ gv_tmp = 0.0_rk
183	!!$ call Vect_Gather(arr_copy, gv_tmp, Msh)
184	!!$ if (ismaster()) &
185	!!$ call Vect_Print(gv_tmp, 'r ::')
186
187	init_norm = Vect_dot_product(r,r)
188	if (init_norm == 0.0) &
189	init_norm = 1.0
190
191	write(stream,'(a,i6,a,e8.2)') 'maxit = ',maxit,', tol = ',tol
192	write(stream,'(a,e22.15)') 'init_norm = ', dsqrt(init_norm)
193
194	iter = 1
195	ratio_norm = 1.0_rk
196	do while((ratio_norm > tol*tol).and.(iter <= maxit))
197	call msolve(m,r,z)
198	!if (present(CoarseMtx_)) then
199	! call preconditioner(z,A,r,CoarseMtx_)
200	! !z=r
201	!else
202	! call preconditioner(z,A,r)
203	!endif
204
205	!call Vect_Print(z,'preconditioner')
206
207	! compute current rho
208	rho_curr = Vect_dot_product(r,z)
209	if (iter == 1) then
210	p = z
211	else
212	beta = rho_curr / rho_prev
213	p = z + beta * p
214	end if
215	!q = SpMtx_pmvm(A, p, Msh)
216	call SpMtx_pmvm(q,A,p,Msh)
217	! compute alpha
218	alpha = rho_curr / Vect_dot_product(p,q)
219	x = x + alpha * p
220	r = r - alpha * q
221	rho_prev = rho_curr
222	! check
223	res_norm = Vect_dot_product(r,r)
224	ratio_norm = dsqrt(res_norm) / dsqrt(init_norm)
225	!!$ ratio_norm = res_norm / init_norm
226	if (ismaster()) &
227	write(stream, '(i5,a,e22.15)') iter,': res_norm=',dsqrt(res_norm)
228	!write(stream, '(i5,a,e22.15,e22.15)') iter,': dsqrt(ratio_norm)=',&
229	!dsqrt(ratio_norm), dsqrt(res_norm)
230	iter = iter + 1
231	end do
232
233	! Deallocate auxiliary data structures
234	! helped to assist with pmvm
235	call pmvmCommStructs_destroy()
236
237	deallocate(b_k)
238	deallocate(z)
239	deallocate(m)
240	deallocate(q)
241	deallocate(p)
242	deallocate(r)
243	if (associated(gv_tmp)) deallocate(gv_tmp)
244	if (associated(arr_copy)) deallocate(arr_copy)
245	end subroutine pcg
246
247	subroutine prec1Level(sol,A,rhs,M,res,A_interf_,CoarseMtx_,Restrict,refactor_)
248	implicit none
249	real(kind=rk),dimension(:),pointer :: sol !< solution
250	type(SpMtx) :: A !< sparse system matrix
251	real(kind=rk),dimension(:),pointer :: rhs !< right hand side
252	type(Mesh),intent(in) :: M !< Mesh
253	real(kind=rk),dimension(:),pointer :: res !< residual vector, allocated
254	!! here for multiplicative Schwarz
255	type(SpMtx),optional :: A_interf_ !< matr@interf.
256	type(SpMtx),optional :: CoarseMtx_ !< Coarse matrix
257	type(SpMtx),optional :: Restrict !< Restriction matrix
258	logical,intent(inout),optional :: refactor_
259
260	type(SpMtx) :: A_tmp
261
262	if(sctls%method==1) then
263	if (refactor_) then!{
264	if (.not.present(CoarseMtx_).or.sctls%input_type==DCTL_INPUT_TYPE_ELEMENTAL.or.numprocs>1) then
265	call Factorise_subdomains(A,A_interf_)
266	else
267	if (.not.present(Restrict)) call DOUG_abort("Restriction matrix needs to be passed along with the coarse matrix!")
268	call Factorise_subdomains(A,AC=CoarseMtx_)
269	end if
270	refactor_=.false.
271	end if
272	call solve_subdomains(sol,A,rhs)
273
274	else if (sctls%method>1) then ! For multiplicative Schwarz method...:
275	if (numprocs>1) call DOUG_abort('multiplicative Schwarz only for numprocs==1 so far',-1)
276	if (.not.present(CoarseMtx_).or.sctls%input_type==DCTL_INPUT_TYPE_ELEMENTAL) then
277	call multiplicative_sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, &
278	A_interf_=A_interf_, &
279	refactor=refactor_) !fine solves
280	else
281	if (.not.present(Restrict)) call DOUG_abort("Restriction matrix needs to be passed along with the coarse matrix!")
282	call multiplicative_sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, &
283	A_interf_=A_interf_,AC=CoarseMtx_, &
284	refactor=refactor_,Restrict=Restrict) !fine solves
285	end if
286	refactor_=.false.
287	endif
288	end subroutine prec1Level
289
290	subroutine prec2Level(prepare,A,sol,rhs,res,CoarseMtx_,Restrict,isFirstIter)
291	use CoarseAllgathers
292
293	logical, intent(in) :: prepare
294	type(SpMtx) :: A !< sparse system matrix
295	real(kind=rk),dimension(:),pointer :: sol !< solution
296	real(kind=rk),dimension(:),pointer :: rhs !< right hand side
297	real(kind=rk),dimension(:),pointer :: res !< residual vector
298	type(SpMtx) :: CoarseMtx_ !< Coarse matrix
299	type(SpMtx) :: Restrict !< Restriction matrix
300	logical :: isFirstIter
301
302	real(kind=rk),dimension(:),pointer,save :: csol,crhs,clrhs,tmpsol2,tmpsol
303	integer :: ol
304
305	ol=max(sctls%overlap,sctls%smoothers)
306
307	if (prepare) then
308	call prec2Level_prepare()
309	else
310	call prec2Level_solve()
311	end if
312
313	contains
314	subroutine prec2Level_exchangeMatrix()
315	integer :: ol
316	logical :: add
317
318	ol=max(sctls%overlap,sctls%smoothers)
319	if (ol==0) then
320	add=.false.
321	else
322	add=.true.
323	endif
324
325	call AllSendCoarseMtx(cdat%LAC,CoarseMtx_,cdat%lg_cfmap,&
326	cdat%ngfc,cdat%nprocs,cdat%send)
327	call AllRecvCoarseMtx(CoarseMtx_,cdat%send,add=add) ! Recieve it
328
329	end subroutine prec2Level_exchangeMatrix
330
331	subroutine prec2Level_prepare()
332	if (isFirstIter) then
333	if (cdat%active) then
334	! First iteration - send matrix
335	call prec2Level_exchangeMatrix()
336	end if
337
338	! after coarse matrix is received allocate coarse vectors
339	if (associated(crhs)) then
340	if (size(crhs)/=CoarseMtx_%ncols) then
341	deallocate(csol)
342	deallocate(crhs)
343	allocate(crhs(CoarseMtx_%ncols))
344	allocate(csol(CoarseMtx_%nrows))
345	endif
346	else
347	allocate(crhs(CoarseMtx_%ncols))
348	allocate(csol(CoarseMtx_%nrows))
349	endif
350	write(stream,*) "Restrict%nrows", Restrict%nrows
351	allocate(clrhs(Restrict%nrows)) ! allocate memory for vector
352	end if
353
354	if (.not.associated(tmpsol)) then
355	!allocate(tmpsol(A%nrows))
356	allocate(tmpsol(size(rhs)))
357	endif
358
359	if (cdat%active) then
360	! Send coarse vector
361	write(stream,*) "Send coarse vector!"
362	call SpMtx_Ax(clrhs,Restrict,rhs,dozero=.true.) ! restrict <RA>
363	if (cdat_vec%active) then
364	call AllSendCoarseVector(clrhs,cdat_vec%nprocs,cdat_vec%cdisps,&
365	cdat_vec%send,useprev=.not.isFirstIter)
366	else
367	call AllSendCoarseVector(clrhs,cdat%nprocs,cdat%cdisps,&
368	cdat%send,useprev=.not.isFirstIter)
369	endif
370	end if ! cdat%active
371	end subroutine prec2Level_prepare
372
373	subroutine prec2Level_solve()
374	if (.not.cdat%active) then ! 1 processor case
375	if (sctls%smoothers==-1) then
376	allocate(tmpsol2(A%nrows))
377	tmpsol2=0.0_rk
378	call exact_sparse_multismoother(tmpsol2,A,rhs)
379	call SpMtx_Ax(crhs,Restrict,tmpsol2,dozero=.true.) ! restriction
380	else
381	if (sctls%method>1.and.sctls%method/=5) then ! multiplicative Schwarz
382	call SpMtx_Ax(crhs,Restrict,res,dozero=.true.) ! restriction
383	else
384	call SpMtx_Ax(crhs,Restrict,rhs,dozero=.true.) ! restriction
385	endif
386	endif
387	end if
388
389	!csol=0.0_rk
390	if (cdat%active) then
391	! Recieve the vector for solve
392	if (cdat_vec%active) then
393	call AllRecvCoarseVector(crhs,cdat_vec%nprocs,&
394	cdat_vec%cdisps,cdat_vec%glg_cfmap,cdat_vec%send)
395	else
396	call AllRecvCoarseVector(crhs,cdat%nprocs,&
397	cdat%cdisps,cdat%glg_cfmap,cdat%send)
398	endif
399	write(stream,*) "Got coarse vector!"
400	!call MPI_BARRIER(MPI_COMM_WORLD,i)
401	end if
402
403	if (isFirstIter) then
404	write (stream,*) &
405	'factorising coarse matrix of size',CoarseMtx_%nrows, &
406	' and nnz:',CoarseMtx_%nnz
407	call free_spmtx_subsolves(CoarseMtx_)
408	allocate(CoarseMtx_%DD%subsolve_ids(1))
409	CoarseMtx_%DD%subsolve_ids=0
410	CoarseMtx_%DD%nsubsolves=1
411	end if
412
413	write (stream,*) "Coarse solve"
414	! Coarse solve
415	call sparse_singlesolve(CoarseMtx_%DD%subsolve_ids(1),csol,crhs,&
416	nfreds=CoarseMtx_%nrows, &
417	nnz=CoarseMtx_%nnz, &
418	indi=CoarseMtx_%indi, &
419	indj=CoarseMtx_%indj, &
420	val=CoarseMtx_%val)
421	if (isFirstIter) then
422	CoarseMtx_%indi=CoarseMtx_%indi+1
423	CoarseMtx_%indj=CoarseMtx_%indj+1
424	end if
425
426	if (cdat_vec%active) then
427	write(stream,*) "clrhs", shape(clrhs), cdat_vec%gl_cfmap
428	call Vect_remap(csol,clrhs,cdat_vec%gl_cfmap,dozero=.true.)
429	write(stream,*) "Restrict%nrows", Restrict%nrows
430	call SpMtx_Ax(tmpsol,Restrict,clrhs,dozero=.true.,transp=.true.)
431	!write(stream,*) "tmpsol", tmpsol
432	elseif (cdat%active) then
433	call Vect_remap(csol,clrhs,cdat%gl_cfmap,dozero=.true.)
434	call SpMtx_Ax(tmpsol,Restrict,clrhs,dozero=.true.,transp=.true.)
435	else
436	if (sctls%smoothers==-1) then
437	call SpMtx_Ax(tmpsol2,Restrict,csol,dozero=.true.,transp=.true.) ! interpolation
438	tmpsol=0.0_rk
439	call exact_sparse_multismoother(tmpsol,A,tmpsol2)
440	else
441	call SpMtx_Ax(tmpsol,Restrict,csol,dozero=.true.,transp=.true.) ! interpolation
442	endif
443	endif
444
445	if (sctls%method==1) then
446	!call Print_Glob_Vect(sol,M,'sol===',chk_endind=M%ninner)
447	sol(1:A%nrows)=sol(1:A%nrows)+tmpsol(1:A%nrows)
448	elseif (sctls%method==3) then ! fully multiplicative Schwarz
449	sol(1:A%nrows)=sol(1:A%nrows)+tmpsol(1:A%nrows)
450	! calculate the residual:
451	call SpMtx_Ax(res,A,sol,dozero=.true.) !
452	res=rhs-res
453	elseif (sctls%method==2.and.ol>0) then ! fully multiplicative Schwarz
454	sol(1:A%nrows)=tmpsol(1:A%nrows)
455	! calculate the residual:
456	call SpMtx_Ax(res,A,sol,dozero=.true.) !
457	res=rhs-res
458	endif
459	if (((ol==0.and.sctls%method==2).or.sctls%method==5).and.sctls%levels>1) then
460	! multiplicative on fine level, additive with coarse level:
461	sol(1:A%nrows)=sol(1:A%nrows)+tmpsol(1:A%nrows)
462	endif
463	end subroutine prec2Level_solve
464
465	end subroutine prec2Level
466
467	!-------------------------------
468	!> Make preconditioner
469	!-------------------------------
470	subroutine preconditioner(sol,A,rhs,M, &
471	A_interf_,CoarseMtx_,Restrict,refactor_,bugtrack_)
472	use CoarseAllgathers
473	use CoarseMtx_mod
474	use Vect_mod
475	implicit none
476	real(kind=rk),dimension(:),pointer :: sol !< solution
477	type(SpMtx) :: A !< sparse system matrix
478	real(kind=rk),dimension(:),pointer :: rhs !< right hand side
479	type(Mesh),intent(in) :: M !< Mesh
480	real(kind=rk),dimension(:),pointer :: res !< residual vector, allocated
481	!! here for multiplicative Schwarz
482	type(SpMtx),optional :: A_interf_ !< matr@interf.
483	type(SpMtx),optional :: CoarseMtx_ !< Coarse matrix
484	type(SpMtx),optional :: Restrict !< Restriction matrix
485	logical,intent(inout),optional :: refactor_
486	logical,optional :: bugtrack_
487	! ----- local: ------
488	integer :: i
489	logical :: add,bugtrack
490	real(kind=rk) :: t1
491	logical :: isFirstIter
492
493	t1 = MPI_WTime()
494
495	if (present(bugtrack_)) then
496	bugtrack=bugtrack_
497	else
498	bugtrack=.false.
499	endif
500	! ----------------------------
501	isFirstIter = .false.
502	if (present(refactor_)) isFirstIter = refactor_
503
504	if (sctls%method==0) then
505	sol=rhs
506	return
507	endif
508
509	sol=0.0_rk
510	if (sctls%method>1) then ! For multiplicative Schwarz method...:
511	allocate(res(size(rhs)))
512	endif
513
514	if (sctls%levels>1) then
515	call prec2Level(.true.,A,sol,rhs,res,CoarseMtx_,Restrict,isFirstIter)
516	end if
517
518	! first level prec
519	call prec1Level(sol,A,rhs,M,res,A_interf_,CoarseMtx_,Restrict,refactor_)
520
521	if (sctls%levels>1) then
522	call prec2Level(.false.,A,sol,rhs,res,CoarseMtx_,Restrict,isFirstIter)
523	end if
524
525	if (sctls%method>1) then ! For multiplicative Schwarz method...:
526	call prec1Level(sol,A,rhs,M,res,A_interf_,CoarseMtx_,Restrict,refactor_)
527	if (associated(res)) deallocate(res)
528	endif
529
530	time_preconditioner = time_preconditioner + MPI_WTime()-t1
531
532	end subroutine preconditioner
533
534	!--------------------------
535	! Solve
536	!--------------------------
537	subroutine msolve (m, r, z)
538	implicit none
539	real(kind=rk), dimension(:), intent(in) :: m, r
540	real(kind=rk), dimension(:), intent(in out) :: z
541	integer :: i
542
543	do i = 1, size(r)
544	z(i) = r(i) !/ m(i)
545	end do
546	end subroutine msolve
547
548	!--------------------------
549	!> Preconditioned conjugent gradient method with eigenvalues
550	!--------------------------
551	subroutine pcg_weigs (A,b,x,Msh,it,cond_num,A_interf_,tol_,maxit_, &
552	x0_,solinf,resvects_,CoarseMtx_,Restrict,refactor_)
553	use CoarseAllgathers
554
555	implicit none
556
557	type(SpMtx),intent(in out) :: A !< System matrix (sparse)
558	float(kind=rk),dimension(:),pointer :: b !< right hand side
559	float(kind=rk),dimension(:),pointer :: x !< Solution
560	type(Mesh),intent(in) :: Msh !< Mesh - aux data for Ax operation
561
562	integer,intent(out) :: it
563	real(kind=rk),intent(out) :: cond_num
564
565	! Optional arguments
566	type(SpMtx),intent(in out),optional :: A_interf_ !< matr@interf.
567	real(kind=rk), intent(in), optional :: tol_ !< Tolerance of the method
568	integer, intent(in), optional :: maxit_ !< Max number of iterations
569	float(kind=rk), dimension(:), intent(in), optional :: x0_ !< Initial guess
570	type(ConvInf), intent(in out), optional :: solinf !< Solution statistics
571	logical, intent(in), optional :: resvects_ !< Fill in the 'resvect' or not
572	type(SpMtx),optional :: CoarseMtx_ !< Coarse matrix
573	type(SpMtx),optional :: Restrict !< Restriction mtx
574	logical,intent(in),optional :: refactor_
575
576	! Local variables
577	logical :: refactor
578
579	real(kind=rk) :: tol ! Tolerance
580	integer :: maxit ! Max number of iterations
581	logical :: resvects=.false.
582
583	real(kind=rk),dimension(:),pointer,save :: r, p, q, m, z, b_k
584	real(kind=rk),dimension(:),pointer,save :: ztmp !TODO remove me
585	real(kind=rk) :: rho_curr, rho_prev, init_norm, res_norm
586	real(kind=rk) :: tmp,ratio_norm
587	integer :: nfreds
588	integer :: ierr
589	real(kind=rk),dimension(:),pointer,save :: dd,ee,alpha,beta
590	!logical,parameter :: bugtrack=.true.
591	logical,parameter :: bugtrack=.false.
592	! profiling
593	real(8) :: time_iterations, t1
594
595	write(stream,'(/a)') 'Preconditioned conjugate gradient:'
596	if (present(refactor_).and.refactor_) then
597	refactor=.true.
598	else
599	refactor=.false.
600	endif
601
602	if (size(b) /= size(x)) &
603	call DOUG_abort('[pcg] : SEVERE : size(b) /= size(x)',-1)
604
605	nfreds=size(b)
606	if (.not.associated(r)) then
607	allocate(r(nfreds))
608	allocate(p(nfreds))
609	allocate(q(nfreds))
610	allocate(m(nfreds))
611	allocate(z(nfreds))
612	allocate(b_k(nfreds))
613	endif
614
615	tol = sctls%solve_tolerance
616	if (present(tol_)) &
617	tol = tol_
618
619	maxit = sctls%solve_maxiters
620	if (maxit==-1) maxit=100000
621	if (present(maxit_)) &
622	maxit = maxit_
623
624	if (present(resvects_).and.(.not.present(solinf))) &
625	call DOUG_abort('[pcg] : SEVERE : "resvects_" must be given along'//&
626	' with "solinf".',-1)
627
628	! Allocate convergence info structure
629	if (present(solinf).and.(.not.present(resvects_))) then
630	call ConvInf_Init(solinf)
631	else if (present(solinf).and.&
632	(present(resvects_).and.(resvects_.eqv.(.true.)))) then
633	call ConvInf_Init(solinf, maxit)
634	resvects = .true.
635	end if
636
637
638	! Initialise auxiliary data structures
639	! to assist with pmvm
640	call pmvmCommStructs_init(A, Msh)
641
642
643	if (bugtrack)call Print_Glob_Vect(x,Msh,'global x===')
644	call SpMtx_pmvm(r,A,x,Msh)
645
646	r = b - r
647	init_norm = Vect_dot_product(r,r)
648	if (init_norm == 0.0) &
649	init_norm = 1.0
650
651
652	write(stream,'(a,i6,a,e8.2)') 'maxit = ',maxit,', tol = ',tol
653	write(stream,'(a,e22.15)') 'init_norm = ', dsqrt(init_norm)
654
655	it = 0
656	ratio_norm = 1.0_rk
657	! arrays for eigenvalue calculation:
658	if (.not.associated(dd)) then
659	allocate(dd(maxit))
660	allocate(ee(maxit))
661	allocate(alpha(maxit))
662	allocate(beta(maxit))
663	endif
664
665	time_iterations = 0.
666	! iterations
667	do while((ratio_norm > tol*tol).and.(it < maxit))
668	it = it + 1
669	t1 = MPI_WTime()
670
671	if (bugtrack)call Print_Glob_Vect(r,Msh,'global r===',chk_endind=Msh%ninner)
672	call preconditioner(sol=z, &
673	A=A, &
674	rhs=r, &
675	M=Msh, &
676	A_interf_=A_interf_, &
677	CoarseMtx_=CoarseMtx_, &
678	Restrict=Restrict, &
679	refactor_=refactor, &
680	bugtrack_=bugtrack)
681
682	refactor=.false.
683	if (bugtrack)call Print_Glob_Vect(z,Msh,'global bef comm z===',chk_endind=Msh%ninner)
684	if (sctls%method/=0) then
685	call Add_common_interf(z,A,Msh)
686	endif
687	!call Print_Glob_Vect(z,Msh,'global aft comm z===',chk_endind=Msh%ninner)
688	!call Print_Glob_Vect(z,Msh,'global z===')
689	! compute current rho
690	rho_curr = Vect_dot_product(r,z)
691
692	if (it == 1) then
693	p = z
694	else
695	beta(it) = rho_curr / rho_prev
696	p = z + beta(it) * p
697	end if
698	call SpMtx_pmvm(q,A, p, Msh)
699	if (bugtrack)call Print_Glob_Vect(q,Msh,'global q===')
700	! compute alpha
701	alpha(it) = rho_curr / Vect_dot_product(p,q)
702	x = x + alpha(it) * p
703	if (bugtrack)call Print_Glob_Vect(x,Msh,'global x===')
704	r = r - alpha(it) * q
705	if (bugtrack)call Print_Glob_Vect(r,Msh,'global r===')
706	rho_prev = rho_curr
707	! check
708	res_norm = Vect_dot_product(r,r)
709	!write (stream,*) "Norm is", res_norm
710	ratio_norm = res_norm / init_norm
711
712	time_iterations = time_iterations + MPI_WTime()-t1
713
714	if (ismaster()) &
715	write(stream, '(i5,a,e22.15)') it,': res_norm=',dsqrt(res_norm)
716	end do
717
718	if(pstream/=0) then
719	write(pstream, "(I0,':pcg iterations:',I0)") myrank, it
720	write(pstream, "(I0,':iterations time:',F0.3)") myrank, time_iterations
721	write(pstream, "(I0,':preconditioner time:',F0.3)") myrank, time_preconditioner
722	end if
723
724	if (ismaster()) then
725	call CalculateEigenvalues(it,dd,ee,alpha,beta,maxit)
726	cond_num=dd(it)/dd(1)
727	write(stream,'(a,i3,a,e10.4)') '#it:',it,' Cond#: ',cond_num
728	endif
729	!deallocate(beta)
730	!deallocate(alpha)
731	!deallocate(ee)
732	!deallocate(dd)
733	! Deallocate auxiliary data structures
734	! helped to assist with pmvm
735	! call pmvmCommStructs_destroy()
736
737	!deallocate(b_k)
738	!deallocate(z)
739	!deallocate(m)
740	!deallocate(q)
741	!deallocate(p)
742	!deallocate(r)
743	end subroutine pcg_weigs
744
745
746	subroutine CalculateEigenvalues(it,dd,ee,alpha,beta,MaxIt)
747	! Finds eigenvalues using Lanczos connection
748	implicit none
749	integer :: it,MaxIt
750	real(kind=rk),dimension(:),pointer :: dd,ee,alpha,beta
751	integer :: i,j,ierr
752
753	dd(1) = 1.0_rk/alpha(1)
754	ee(1) = 0.0_rk
755	do i=2,it
756	if (alpha(i)==0.or.beta(i)<0) then
757	print *,'Unable to compute eigenvalue number',i,' and '
758	do j=1,i
759	print *,j,' alpha:',alpha(j),' beta:',beta
760	enddo
761	return
762	else
763	dd(i) = 1.0_rk/alpha(i) + beta(i)/alpha(i-1)
764	ee(i) = -dsqrt(beta(i))/alpha(i-1)
765	endif
766	enddo
767	!call tql1(it,dd,ee,MaxIt,ierr)
768	call tql1(it,dd,ee,ierr)
769	if (ierr > 0) then
770	print *,'Cancelled after EigMaxIt while calculating eig',ierr
771	endif
772	end subroutine CalculateEigenvalues
773
774	subroutine tql1(n,d,e,ierr)
775	implicit none
776	integer :: i,j,l,m,n,ii,l1,l2,mml,ierr
777	real(kind=rk) :: d(n),e(n)
778	real(kind=rk) :: c,c2,c3,dl1,el1,f,g,h,p,r,s,s2,tst1,tst2 !,pythag
779	!
780	! this subroutine is a translation of the algol procedure tql1,
781	! num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and
782	! wilkinson.
783	! handbook for auto. comp., vol.ii-linear algebra, 227-240(1971).
784	!
785	! this subroutine finds the eigenvalues of a symmetric
786	! tridiagonal matrix by the ql method.
787	!
788	! on input
789	! n is the order of the matrix.
790	! d contains the diagonal elements of the input matrix.
791	! e contains the subdiagonal elements of the input matrix
792	! in its last n-1 positions. e(1) is arbitrary.
793	!
794	! on output
795	! d contains the eigenvalues in ascending order. if an
796	! error exit is made, the eigenvalues are correct and
797	! ordered for indices 1,2,...ierr-1, but may not be
798	! the smallest eigenvalues.
799	! e has been destroyed.
800	! ierr is set to
801	! zero for normal return,
802	! j if the j-th eigenvalue has not been
803	! determined after 30 iterations.
804	! calls pythag for sqrt(aa + bb) .
805	!
806	! questions and comments should be directed to burton s. garbow,
807	! mathematics and computer science div, argonne national laboratory
808	!
809	! this version dated august 1983.
810	!
811	! ------------------------------------------------------------------
812	!
813	ierr = 0
814	if (n .eq. 1) go to 1001
815	!
816	do 100 i = 2, n
817	100 e(i-1) = e(i)
818	!
819	f = 0.0e0
820	tst1 = 0.0e0
821	e(n) = 0.0e0
822
823	do 290 l = 1, n
824	j = 0
825	h = abs(d(l)) + abs(e(l))
826	if (tst1 .lt. h) tst1 = h
827	! .......... look for small sub-diagonal element ..........
828	do 110 m = l, n
829	tst2 = tst1 + abs(e(m))
830	if (tst2 .eq. tst1) go to 120
831	! .......... e(n) is always zero, so there is no exit
832	! throu gh the bottom of the loop ..........
833	110 continue
834
835	120 if (m .eq. l) go to 210
836	130 if (j .eq. 30 ) go to 1000
837	j = j + 1
838	! .......... form shift ..........
839	l1 = l + 1
840	l2 = l1 + 1
841	g = d(l)
842	p = (d(l1) - g) / (2.0e0 * e(l))
843	r = pythag(p,1.0e0_rk)
844	d(l) = e(l) / (p + sign(r,p))
845	d(l1) = e(l) * (p + sign(r,p))
846	dl1 = d(l1)
847	h = g - d(l)
848	if (l2 .gt. n) go to 145
849
850	do 140 i = l2, n
851	140 d(i) = d(i) - h
852
853	145 f = f + h
854	! .......... ql transformation ..........
855	p = d(m)
856	c = 1.0e0
857	c2 = c
858	el1 = e(l1)
859	s = 0.0e0
860	mml = m - l
861	! .......... for i=m-1 step -1 until l do -- ..........
862	do 200 ii = 1, mml
863	c3 = c2
864	c2 = c
865	s2 = s
866	i = m - ii
867	g = c * e(i)
868	h = c * p
869	r = pythag(p,e(i))
870	e(i+1) = s * r
871	s = e(i) / r
872	c = p / r
873	p = c * d(i) - s * g
874	d(i+1) = h + s * (c * g + s * d(i))
875	200 continue
876
877	p = -s * s2 * c3 * el1 * e(l) / dl1
878	e(l) = s * p
879	d(l) = c * p
880	tst2 = tst1 + abs(e(l))
881	if (tst2 .gt. tst1) go to 130
882	210 p = d(l) + f
883	! .......... order eigenvalues ..........
884	if (l .eq. 1) go to 250
885	! .......... for i=l step -1 until 2 do -- ..........
886	do 230 ii = 2, l
887	i = l + 2 - ii
888	if (p .ge. d(i-1)) go to 270
889	d(i) = d(i-1)
890	230 continue
891
892	250 i = 1
893	270 d(i) = p
894	290 continue
895
896	go to 1001
897	! .......... set error -- no convergence to an
898	! eigenvalue after 30 iterations ..........
899	1000 ierr = l
900	1001 return
901	end subroutine tql1
902
903	function pythag(a,b) result(pyt)
904	! finds dsqrt(a2+b2) without overflow or destructive underflow
905	implicit none
906	real(kind=rk) :: a,b,pyt
907	real(kind=rk) :: p,r,s,t,u
908
909	p = dmax1(dabs(a),dabs(b))
910	if (p .eq. 0.0d0) go to 20
911	r = (dmin1(dabs(a),dabs(b))/p)**2
912	10 continue
913	t = 4.0d0 + r
914	if (t .eq. 4.0d0) go to 20
915	s = r/t
916	u = 1.0d0 + 2.0d0*s
917	p = u*p
918	r = (s/u)*2 r
919	go to 10
920	20 pyt = p
921	return
922	end function pythag
923
924
925	end module pcg_mod
926
927

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

Download in other formats:

Original Format