Context Navigation

pcg.F90 @ 82697a3

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

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

Download in other formats:

Original Format