Context Navigation

pcg.F90 @ 7efae20

Revision 7efae20, 35.4 KB checked in by chris <chris@…>, 6 years ago (diff)

more doxygen comments

git-svn-id: svn://dougdevel.org/doug/trunk@212 146eac83-2b0f-0410-9363-9ebcc3d868b9

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	!-------------------------------
246	!> Make preconditioner
247	!-------------------------------
248	subroutine preconditioner(sol,A,rhs,M, &
249	A_interf_,CoarseMtx_,Restrict,refactor_,bugtrack_)
250	use subsolvers
251	use CoarseMtx_mod
252	use CoarseAllgathers
253	use Vect_mod
254	implicit none
255	real(kind=rk),dimension(:),pointer :: sol !< solution
256	type(SpMtx) :: A !< sparse system matrix
257	real(kind=rk),dimension(:),pointer :: rhs !< right hand side
258	type(Mesh),intent(in) :: M !< Mesh
259	real(kind=rk),dimension(:),pointer :: res !< residual vector, allocated
260	!! here for multiplicative Schwarz
261	type(SpMtx),optional :: A_interf_ !< matr@interf.
262	type(SpMtx),optional :: CoarseMtx_ !< Coarse matrix
263	type(SpMtx),optional :: Restrict !< Restriction matrix
264	logical,intent(inout),optional :: refactor_
265	logical,optional :: bugtrack_
266	! ----- local: ------
267	real(kind=rk),dimension(:),pointer,save :: csol,crhs,tmpsol,tmpsol2,clrhs
268	type(SpMtx) :: A_tmp
269	integer :: i,ol
270	logical :: add,bugtrack
271	real(kind=rk) :: t1
272
273	t1 = MPI_WTime()
274
275	if (present(bugtrack_)) then
276	bugtrack=bugtrack_
277	else
278	bugtrack=.false.
279	endif
280	! ----------------------------
281	if (sctls%method==0) then
282	sol=rhs
283	return
284	endif
285	sol=0.0_rk
286	if (sctls%method>1) then ! For multiplicative Schwarz method...:
287	allocate(res(size(rhs)))
288	endif
289	if (present(CoarseMtx_)) then !{
290	if (.not.present(Restrict)) call DOUG_abort("Restriction matrix needs to be passed along with the coarse matrix!")
291	if (cdat%active) then
292	! Not the first iteration, so send the coarse vector here
293	if (.not.(present(refactor_).and.refactor_)) then
294	call SpMtx_Ax(clrhs,Restrict,rhs,dozero=.true.) ! restrict <RA>
295	! And send
296	if (bugtrack)write(stream,*) "local coarse RHS is:",clrhs
297	if (cdat_vec%active) then
298	call AllSendCoarseVector(clrhs,cdat_vec%nprocs,cdat_vec%cdisps,&
299	cdat_vec%send,useprev=.true.)
300	else
301	call AllSendCoarseVector(clrhs,cdat%nprocs,cdat%cdisps,&
302	cdat%send,useprev=.true.)
303	endif
304	else ! First iteration - send matrix
305	call AllSendCoarseMtx(cdat%LAC,CoarseMtx_,cdat%lg_cfmap,&
306	cdat%ngfc,cdat%nprocs,cdat%send)
307	endif
308	endif
309
310	ol=max(sctls%overlap,sctls%smoothers)
311	if (sctls%input_type==DCTL_INPUT_TYPE_ELEMENTAL.or.numprocs>1) then
312	call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, &
313	A_interf_=A_interf_, &
314	refactor=refactor_) !fine solves
315	else
316	call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, &
317	A_interf_=A_interf_,AC=CoarseMtx_, &
318	refactor=refactor_,Restrict=Restrict) !fine solves
319	if (bugtrack)call Print_Glob_Vect(sol,M,'global sol===')
320	endif
321
322	if (sctls%levels>1) then
323	if (ol==0) then
324	add=.false.
325	else
326	add=.true.
327	endif
328	if (present(refactor_).and.refactor_) then !{ setup coarse solve:
329	! We have a matrix send pending
330	if (cdat%active) then
331	!write(stream,*) "Waiting the Matrix!"
332	call AllRecvCoarseMtx(CoarseMtx_,cdat%send,add=add) ! Recieve it
333	!write(stream,*) "Got Matrix!"
334	allocate(clrhs(Restrict%nrows)) ! allocate memory for vector
335	call SpMtx_Ax(clrhs,Restrict,rhs,dozero=.true.) ! restrict <RA>
336	! And send
337	if (bugtrack)write(stream,*) "(f) local Coarse RHS is:",clrhs
338	if (cdat_vec%active) then
339	call AllSendCoarseVector(clrhs,cdat_vec%nprocs,&
340	cdat_vec%cdisps,cdat_vec%send)
341	else
342	call AllSendCoarseVector(clrhs,cdat%nprocs,&
343	cdat%cdisps,cdat%send)
344	endif
345	endif
346	if (associated(crhs)) then
347	if (size(crhs)/=CoarseMtx_%ncols) then
348	deallocate(csol)
349	deallocate(crhs)
350	allocate(crhs(CoarseMtx_%ncols))
351	allocate(csol(CoarseMtx_%nrows))
352	endif
353	else
354	allocate(crhs(CoarseMtx_%ncols))
355	allocate(csol(CoarseMtx_%nrows))
356	endif
357	if (.not.associated(tmpsol)) then
358	!allocate(tmpsol(A%nrows))
359	allocate(tmpsol(size(rhs)))
360	endif
361
362	if (.not.cdat%active) then
363	if (sctls%smoothers==-1) then
364	allocate(tmpsol2(A%nrows))
365	tmpsol2=0.0_rk
366	call exact_sparse_multismoother(tmpsol2,A,rhs)
367	call SpMtx_Ax(crhs,Restrict,tmpsol2,dozero=.true.) ! restriction
368	else
369	if (sctls%method>1.and.sctls%method/=5) then ! multiplicative Schwarz
370	call SpMtx_Ax(crhs,Restrict,res,dozero=.true.) ! restriction
371	else
372	call SpMtx_Ax(crhs,Restrict,rhs,dozero=.true.) ! restriction
373	endif
374	endif
375	endif
376	!csol=0.0_rk
377	write (stream,*) &
378	'factorising coarse matrix of size',CoarseMtx_%nrows, &
379	' and nnz:',CoarseMtx_%nnz
380	call free_spmtx_subsolves(CoarseMtx_)
381	allocate(CoarseMtx_%subsolve_ids(1))
382	CoarseMtx_%subsolve_ids=0
383	CoarseMtx_%nsubsolves=1
384
385	!! TODO: force it to do it AFTER factorization, but before solve
386	!if (present(cdat)) then
387	if (cdat%active) then
388	! Factorise the matrix
389	if (sctls%verbose>2) then
390	write(stream,*)'Global coarse matrix is:---------'
391	do i=1,CoarseMtx_%nnz
392	write(stream,*) CoarseMtx_%indi(i),&
393	CoarseMtx_%indj(i),CoarseMtx_%val(i)
394	enddo
395	write(stream,*)'---------------------------------'
396	endif
397	call factorise(CoarseMtx_%subsolve_ids(1), &
398	nfreds=CoarseMtx_%nrows, &
399	nnz=CoarseMtx_%nnz, &
400	indi=CoarseMtx_%indi, &
401	indj=CoarseMtx_%indj, &
402	val=CoarseMtx_%val)
403	! Recieve the vector for solve
404	if (cdat_vec%active) then
405	call AllRecvCoarseVector(crhs,cdat_vec%nprocs,&
406	cdat_vec%cdisps,cdat_vec%glg_cfmap,cdat_vec%send)
407	else
408	call AllRecvCoarseVector(crhs,cdat%nprocs,&
409	cdat%cdisps,cdat%glg_cfmap,cdat%send)
410	endif
411	!write(stream,*) "Got coarse vector!"
412	!call MPI_BARRIER(MPI_COMM_WORLD,i)
413	endif
414	if (bugtrack) then
415	if (isnan(sum(crhs))) then
416	write(stream,*) "(f) global Coarse RHS is:",crhs
417	call DOUG_abort('Got NaN receive fromAllRecvCoarseVector...',838465)
418	endif
419	!else
420	! write(stream,*) "(f) global Coarse RHS is:",crhs
421	endif
422
423
424	call sparse_singlesolve(CoarseMtx_%subsolve_ids(1),csol,crhs, &
425	nfreds=CoarseMtx_%nrows, &
426	nnz=CoarseMtx_%nnz, &
427	indi=CoarseMtx_%indi, &
428	indj=CoarseMtx_%indj, &
429	val=CoarseMtx_%val)
430	write (stream,*) 'coarse factorisation done!',CoarseMtx_%subsolve_ids(1)
431	if (bugtrack)write(stream,*) "(f) Coarse SOL is:",csol
432
433	if (cdat_vec%active) then
434	CoarseMtx_%indi=CoarseMtx_%indi+1
435	CoarseMtx_%indj=CoarseMtx_%indj+1
436	call Vect_remap(csol,clrhs,cdat_vec%gl_cfmap,dozero=.true.)
437	call SpMtx_Ax(tmpsol,Restrict,clrhs,dozero=.true.,transp=.true.)
438	elseif (cdat%active) then
439	CoarseMtx_%indi=CoarseMtx_%indi+1
440	CoarseMtx_%indj=CoarseMtx_%indj+1
441	call Vect_remap(csol,clrhs,cdat%gl_cfmap,dozero=.true.)
442	call SpMtx_Ax(tmpsol,Restrict,clrhs,dozero=.true.,transp=.true.)
443	else
444	if (sctls%smoothers==-1) then
445	call SpMtx_Ax(tmpsol2,Restrict,csol,dozero=.true.,transp=.true.) ! interpolation
446	tmpsol=0.0_rk
447	call exact_sparse_multismoother(tmpsol,A,tmpsol2)
448	else
449	call SpMtx_Ax(tmpsol,Restrict,csol,dozero=.true.,transp=.true.) ! interpolation
450	endif
451	endif
452	if (bugtrack)call Print_Glob_Vect(tmpsol,M,'tmpsol===',chk_endind=M%ninner)
453	else !}{ apply coarse solve:
454	if (cdat_vec%active) then
455	call AllRecvCoarseVector(crhs,cdat_vec%nprocs,cdat_vec%cdisps,&
456	cdat_vec%glg_cfmap,cdat_vec%send)
457	elseif (cdat%active) then
458	call AllRecvCoarseVector(crhs,cdat%nprocs,cdat%cdisps,&
459	cdat%glg_cfmap,cdat%send)
460	else
461	if (sctls%smoothers==-1) then
462	tmpsol2=0.0_rk
463	call exact_sparse_multismoother(tmpsol2,A,rhs)
464	call SpMtx_Ax(crhs,Restrict,tmpsol2,dozero=.true.) ! restriction
465	else
466	call SpMtx_Ax(crhs,Restrict,rhs,dozero=.true.) ! restriction
467	endif
468	endif
469	call sparse_singlesolve(CoarseMtx_%subsolve_ids(1),csol,crhs, &
470	nfreds=CoarseMtx_%nrows)
471	if (cdat_vec%active) then
472	call Vect_remap(csol,clrhs,cdat_vec%gl_cfmap,dozero=.true.)
473	call SpMtx_Ax(tmpsol,Restrict,clrhs,dozero=.true.,transp=.true.)
474	elseif (cdat%active) then
475	call Vect_remap(csol,clrhs,cdat%gl_cfmap,dozero=.true.)
476	call SpMtx_Ax(tmpsol,Restrict,clrhs,dozero=.true.,transp=.true.)
477	else
478	if (sctls%smoothers==-1) then
479	call SpMtx_Ax(tmpsol2,Restrict,csol,dozero=.true.,transp=.true.) ! interpolation
480	tmpsol=0.0_rk
481	call exact_sparse_multismoother(tmpsol,A,tmpsol2)
482	else
483	call SpMtx_Ax(tmpsol,Restrict,csol,dozero=.true.,transp=.true.) ! interpolation
484	! call Vect_print(tmpsol,"csol")
485	endif
486	endif
487	endif !}
488	if (sctls%method==1) then
489	!call Print_Glob_Vect(sol,M,'sol===',chk_endind=M%ninner)
490	sol(1:A%nrows)=sol(1:A%nrows)+tmpsol(1:A%nrows)
491	elseif (sctls%method==3) then ! fully multiplicative Schwarz
492	sol(1:A%nrows)=sol(1:A%nrows)+tmpsol(1:A%nrows)
493	! calculate the residual:
494	call SpMtx_Ax(res,A,sol,dozero=.true.) !
495	res=rhs-res
496	elseif (sctls%method==2.and.ol>0) then ! fully multiplicative Schwarz
497	sol(1:A%nrows)=tmpsol(1:A%nrows)
498	! calculate the residual:
499	call SpMtx_Ax(res,A,sol,dozero=.true.) !
500	res=rhs-res
501	endif
502	endif
503	if (sctls%method>1) then ! For multiplicative Schwarz method...:
504	refactor_=.false.
505	if (sctls%input_type==DCTL_INPUT_TYPE_ELEMENTAL) then
506	call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, &
507	A_interf_=A_interf_, &
508	refactor=refactor_) !fine solves
509	else
510	call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, &
511	A_interf_=A_interf_,AC=CoarseMtx_, &
512	refactor=refactor_,Restrict=Restrict) !fine solves
513	endif
514	endif
515	if (((ol==0.and.sctls%method==2).or.sctls%method==5).and.sctls%levels>1) then
516	! multiplicative on fine level, additive with coarse level:
517	sol(1:A%nrows)=sol(1:A%nrows)+tmpsol(1:A%nrows)
518	endif
519	else !}{
520	if (refactor_.and.present(A_interf_)) then!{
521	if (sctls%verbose>9) then
522	!call SpMtx_printMat(A)
523	call SpMtx_printRaw(A)
524	!call SpMtx_printMat(A_interf_)
525	call SpMtx_printRaw(A_interf_)
526	!stop
527	endif
528	if (sctls%input_type==DCTL_INPUT_TYPE_ASSEMBLED) then
529	! we need the fully sorted entries but A has block struct...
530	! keep the original values...
531	if (associated(A%M_bound)) then
532	A_tmp=SpMtx_newInit(nnz=A%nnz,nblocks=A%nblocks,&
533	nrows=A%nrows,ncols=A%ncols,&
534	indi=A%indi,indj=A%indj,val=A%val,&
535	arrange_type=A%arrange_type,M_bound=A%M_bound)
536	else
537	A_tmp=SpMtx_newInit(nnz=A%nnz,nblocks=A%nblocks,&
538	nrows=A%nrows,ncols=A%ncols,&
539	indi=A%indi,indj=A%indj,val=A%val,&
540	arrange_type=A%arrange_type)
541	endif
542	A%arrange_type=D_SpMtx_ARRNG_NO
543	!A%nrows=max(A%nrows,A_interf_%nrows)
544	!A%ncols=max(A%nrows,A_interf_%ncols)
545	!if (associated(A%M_bound)) then
546	! deallocate(A%M_bound)
547	!endif
548	call SpMtx_arrange(A,D_SpMtx_ARRNG_ROWS,sort=.true.)
549	call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, &
550	A_interf_=A_interf_, &
551	refactor=refactor_) !fine solves
552	if (sctls%method>1) then ! For multiplicative Schwarz method...:
553	refactor_=.false.
554	call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, &
555	A_interf_=A_interf_, &
556	refactor=refactor_) !fine solves
557	endif
558	! put the original structure and orders back:
559	A%indi=A_tmp%indi
560	A%indj=A_tmp%indj
561	A%val=A_tmp%val
562	if (associated(A_tmp%M_bound)) A%M_bound=A_tmp%M_bound
563	A%nrows=A_tmp%nrows
564	A%ncols=A_tmp%ncols
565	A%arrange_type=A_tmp%arrange_type
566	call SpMtx_Destroy(A_tmp)
567	else
568	call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, &
569	A_interf_=A_interf_, &
570	refactor=refactor_) !fine solves
571	if (sctls%method>1) then ! For multiplicative Schwarz method...:
572	refactor_=.false.
573	call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, &
574	A_interf_=A_interf_, &
575	refactor=refactor_) !fine solves
576	endif
577	endif
578	elseif (refactor_.and.sctls%input_type==DCTL_INPUT_TYPE_ASSEMBLED) then!}{
579	if (associated(A%M_bound)) then
580	A_tmp=SpMtx_newInit(nnz=size(A%indi),nblocks=A%nblocks,&
581	nrows=A%nrows,ncols=A%ncols,&
582	indi=A%indi,indj=A%indj,val=A%val,&
583	arrange_type=A%arrange_type,M_bound=A%M_bound)
584	else
585	A_tmp=SpMtx_newInit(nnz=size(A%indi),nblocks=A%nblocks,&
586	nrows=A%nrows,ncols=A%ncols,&
587	indi=A%indi,indj=A%indj,val=A%val,&
588	arrange_type=A%arrange_type)
589	endif
590	A_tmp%nnz=A%nnz
591	A%arrange_type=D_SpMtx_ARRNG_NO
592	call SpMtx_arrange(A,D_SpMtx_ARRNG_ROWS,sort=.true.)
593	call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, &
594	refactor=refactor_) !fine solves
595	if (sctls%method>1) then ! For multiplicative Schwarz method...:
596	refactor_=.false.
597	call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, &
598	refactor=refactor_) !fine solves
599	endif
600	! put the original structure and orders back:
601	A%indi=A_tmp%indi
602	A%indj=A_tmp%indj
603	A%val=A_tmp%val
604	if (associated(A_tmp%M_bound)) A%M_bound=A_tmp%M_bound
605	A%nrows=A_tmp%nrows
606	A%ncols=A_tmp%ncols
607	A%arrange_type=A_tmp%arrange_type
608	call SpMtx_Destroy(A_tmp)
609	else!}{
610	call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, &
611	A_interf_=A_interf_, &
612	refactor=refactor_) !fine solves
613	if (sctls%method>1) then ! For multiplicative Schwarz method...:
614	refactor_=.false.
615	call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, &
616	A_interf_=A_interf_, &
617	refactor=refactor_) !fine solves
618	endif
619	endif!}
620	endif !}
621	if (sctls%method>1) then ! For multiplicative Schwarz method...:
622	if (associated(res)) deallocate(res)
623	endif
624
625	time_preconditioner = time_preconditioner + MPI_WTime()-t1
626
627	end subroutine preconditioner
628
629	!--------------------------
630	! Solve
631	!--------------------------
632	subroutine msolve (m, r, z)
633	implicit none
634	real(kind=rk), dimension(:), intent(in) :: m, r
635	real(kind=rk), dimension(:), intent(in out) :: z
636	integer :: i
637
638	do i = 1, size(r)
639	z(i) = r(i) !/ m(i)
640	end do
641	end subroutine msolve
642
643	!--------------------------
644	!> Preconditioned conjugent gradient method with eigenvalues
645	!--------------------------
646	subroutine pcg_weigs (A,b,x,Msh,it,cond_num,A_interf_,tol_,maxit_, &
647	x0_,solinf,resvects_,CoarseMtx_,Restrict,refactor_)
648	use CoarseAllgathers
649
650	implicit none
651
652	type(SpMtx),intent(in out) :: A !< System matrix (sparse)
653	float(kind=rk),dimension(:),pointer :: b !< right hand side
654	float(kind=rk),dimension(:),pointer :: x !< Solution
655	type(Mesh),intent(in) :: Msh !< Mesh - aux data for Ax operation
656
657	integer,intent(out) :: it
658	real(kind=rk),intent(out) :: cond_num
659
660	! Optional arguments
661	type(SpMtx),intent(in out),optional :: A_interf_ !< matr@interf.
662	real(kind=rk), intent(in), optional :: tol_ !< Tolerance of the method
663	integer, intent(in), optional :: maxit_ !< Max number of iterations
664	float(kind=rk), dimension(:), intent(in), optional :: x0_ !< Initial guess
665	type(ConvInf), intent(in out), optional :: solinf !< Solution statistics
666	logical, intent(in), optional :: resvects_ !< Fill in the 'resvect' or not
667	type(SpMtx),optional :: CoarseMtx_ !< Coarse matrix
668	type(SpMtx),optional :: Restrict !< Restriction mtx
669	logical,intent(in),optional :: refactor_
670
671	! Local variables
672	logical :: refactor
673
674	real(kind=rk) :: tol ! Tolerance
675	integer :: maxit ! Max number of iterations
676	logical :: resvects=.false.
677
678	real(kind=rk),dimension(:),pointer,save :: r, p, q, m, z, b_k
679	real(kind=rk),dimension(:),pointer,save :: ztmp !TODO remove me
680	real(kind=rk) :: rho_curr, rho_prev, init_norm, res_norm
681	real(kind=rk) :: tmp,ratio_norm
682	integer :: nfreds
683	integer :: ierr
684	real(kind=rk),dimension(:),pointer,save :: dd,ee,alpha,beta
685	!logical,parameter :: bugtrack=.true.
686	logical,parameter :: bugtrack=.false.
687
688	write(stream,'(/a)') 'Preconditioned conjugate gradient:'
689	if (present(refactor_).and.refactor_) then
690	refactor=.true.
691	else
692	refactor=.false.
693	endif
694
695	if (size(b) /= size(x)) &
696	call DOUG_abort('[pcg] : SEVERE : size(b) /= size(x)',-1)
697
698	nfreds=size(b)
699	if (.not.associated(r)) then
700	allocate(r(nfreds))
701	allocate(p(nfreds))
702	allocate(q(nfreds))
703	allocate(m(nfreds))
704	allocate(z(nfreds))
705	allocate(b_k(nfreds))
706	endif
707
708	tol = sctls%solve_tolerance
709	if (present(tol_)) &
710	tol = tol_
711
712	maxit = sctls%solve_maxiters
713	if (maxit==-1) maxit=100000
714	if (present(maxit_)) &
715	maxit = maxit_
716
717	if (present(resvects_).and.(.not.present(solinf))) &
718	call DOUG_abort('[pcg] : SEVERE : "resvects_" must be given along'//&
719	' with "solinf".',-1)
720
721	! Allocate convergence info structure
722	if (present(solinf).and.(.not.present(resvects_))) then
723	call ConvInf_Init(solinf)
724	else if (present(solinf).and.&
725	(present(resvects_).and.(resvects_.eqv.(.true.)))) then
726	call ConvInf_Init(solinf, maxit)
727	resvects = .true.
728	end if
729
730
731	! Initialise auxiliary data structures
732	! to assist with pmvm
733	call pmvmCommStructs_init(A, Msh)
734
735
736	if (bugtrack)call Print_Glob_Vect(x,Msh,'global x===')
737	call SpMtx_pmvm(r,A,x,Msh)
738
739	r = b - r
740	init_norm = Vect_dot_product(r,r)
741	if (init_norm == 0.0) &
742	init_norm = 1.0
743
744
745	write(stream,'(a,i6,a,e8.2)') 'maxit = ',maxit,', tol = ',tol
746	write(stream,'(a,e22.15)') 'init_norm = ', dsqrt(init_norm)
747
748	it = 0
749	ratio_norm = 1.0_rk
750	! arrays for eigenvalue calculation:
751	if (.not.associated(dd)) then
752	allocate(dd(maxit))
753	allocate(ee(maxit))
754	allocate(alpha(maxit))
755	allocate(beta(maxit))
756	endif
757	do while((ratio_norm > tol*tol).and.(it < maxit))
758	it = it + 1
759
760	if (bugtrack)call Print_Glob_Vect(r,Msh,'global r===',chk_endind=Msh%ninner)
761	write(6,*) '############# CALLING PRECOND #############' !XXX
762	call preconditioner(sol=z, &
763	A=A, &
764	rhs=r, &
765	M=Msh, &
766	A_interf_=A_interf_, &
767	CoarseMtx_=CoarseMtx_, &
768	Restrict=Restrict, &
769	refactor_=refactor, &
770	bugtrack_=bugtrack)
771	refactor=.false.
772	if (bugtrack)call Print_Glob_Vect(z,Msh,'global bef comm z===',chk_endind=Msh%ninner)
773	if (sctls%method/=0) then
774	call Add_common_interf(z,A,Msh)
775	endif
776	!call Print_Glob_Vect(z,Msh,'global aft comm z===',chk_endind=Msh%ninner)
777	!call Print_Glob_Vect(z,Msh,'global z===')
778	! compute current rho
779	rho_curr = Vect_dot_product(r,z)
780
781	if (it == 1) then
782	p = z
783	else
784	beta(it) = rho_curr / rho_prev
785	p = z + beta(it) * p
786	end if
787	call SpMtx_pmvm(q,A, p, Msh)
788	if (bugtrack)call Print_Glob_Vect(q,Msh,'global q===')
789	! compute alpha
790	alpha(it) = rho_curr / Vect_dot_product(p,q)
791	x = x + alpha(it) * p
792	if (bugtrack)call Print_Glob_Vect(x,Msh,'global x===')
793	r = r - alpha(it) * q
794	if (bugtrack)call Print_Glob_Vect(r,Msh,'global r===')
795	rho_prev = rho_curr
796	! check
797	res_norm = Vect_dot_product(r,r)
798	!write (stream,*) "Norm is", res_norm
799	ratio_norm = res_norm / init_norm
800	if (ismaster()) &
801	write(stream, '(i5,a,e22.15)') it,': res_norm=',dsqrt(res_norm)
802	end do
803
804	if(pstream/=0) then
805	write(pstream, "(I0,':pcg iterations:',I0)") myrank, it
806	write(pstream, "(I0,':preconditioner time:',F0.3)") myrank, time_preconditioner
807	end if
808
809	if (ismaster()) then
810	call CalculateEigenvalues(it,dd,ee,alpha,beta,maxit)
811	cond_num=dd(it)/dd(1)
812	write(stream,'(a,i3,a,e10.4)') '#it:',it,' Cond#: ',cond_num
813	endif
814	!deallocate(beta)
815	!deallocate(alpha)
816	!deallocate(ee)
817	!deallocate(dd)
818	! Deallocate auxiliary data structures
819	! helped to assist with pmvm
820	! call pmvmCommStructs_destroy()
821
822	!deallocate(b_k)
823	!deallocate(z)
824	!deallocate(m)
825	!deallocate(q)
826	!deallocate(p)
827	!deallocate(r)
828	end subroutine pcg_weigs
829
830
831	subroutine CalculateEigenvalues(it,dd,ee,alpha,beta,MaxIt)
832	! Finds eigenvalues using Lanczos connection
833	implicit none
834	integer :: it,MaxIt
835	real(kind=rk),dimension(:),pointer :: dd,ee,alpha,beta
836	integer :: i,j,ierr
837
838	dd(1) = 1.0_rk/alpha(1)
839	ee(1) = 0.0_rk
840	do i=2,it
841	if (alpha(i)==0.or.beta(i)<0) then
842	print *,'Unable to compute eigenvalue number',i,' and '
843	do j=1,i
844	print *,j,' alpha:',alpha(j),' beta:',beta
845	enddo
846	return
847	else
848	dd(i) = 1.0_rk/alpha(i) + beta(i)/alpha(i-1)
849	ee(i) = -dsqrt(beta(i))/alpha(i-1)
850	endif
851	enddo
852	!call tql1(it,dd,ee,MaxIt,ierr)
853	call tql1(it,dd,ee,ierr)
854	if (ierr > 0) then
855	print *,'Cancelled after EigMaxIt while calculating eig',ierr
856	endif
857	end subroutine CalculateEigenvalues
858
859	subroutine tql1(n,d,e,ierr)
860	implicit none
861	integer :: i,j,l,m,n,ii,l1,l2,mml,ierr
862	real(kind=rk) :: d(n),e(n)
863	real(kind=rk) :: c,c2,c3,dl1,el1,f,g,h,p,r,s,s2,tst1,tst2 !,pythag
864	!
865	! this subroutine is a translation of the algol procedure tql1,
866	! num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and
867	! wilkinson.
868	! handbook for auto. comp., vol.ii-linear algebra, 227-240(1971).
869	!
870	! this subroutine finds the eigenvalues of a symmetric
871	! tridiagonal matrix by the ql method.
872	!
873	! on input
874	! n is the order of the matrix.
875	! d contains the diagonal elements of the input matrix.
876	! e contains the subdiagonal elements of the input matrix
877	! in its last n-1 positions. e(1) is arbitrary.
878	!
879	! on output
880	! d contains the eigenvalues in ascending order. if an
881	! error exit is made, the eigenvalues are correct and
882	! ordered for indices 1,2,...ierr-1, but may not be
883	! the smallest eigenvalues.
884	! e has been destroyed.
885	! ierr is set to
886	! zero for normal return,
887	! j if the j-th eigenvalue has not been
888	! determined after 30 iterations.
889	! calls pythag for sqrt(aa + bb) .
890	!
891	! questions and comments should be directed to burton s. garbow,
892	! mathematics and computer science div, argonne national laboratory
893	!
894	! this version dated august 1983.
895	!
896	! ------------------------------------------------------------------
897	!
898	ierr = 0
899	if (n .eq. 1) go to 1001
900	!
901	do 100 i = 2, n
902	100 e(i-1) = e(i)
903	!
904	f = 0.0e0
905	tst1 = 0.0e0
906	e(n) = 0.0e0
907
908	do 290 l = 1, n
909	j = 0
910	h = abs(d(l)) + abs(e(l))
911	if (tst1 .lt. h) tst1 = h
912	! .......... look for small sub-diagonal element ..........
913	do 110 m = l, n
914	tst2 = tst1 + abs(e(m))
915	if (tst2 .eq. tst1) go to 120
916	! .......... e(n) is always zero, so there is no exit
917	! throu gh the bottom of the loop ..........
918	110 continue
919
920	120 if (m .eq. l) go to 210
921	130 if (j .eq. 30 ) go to 1000
922	j = j + 1
923	! .......... form shift ..........
924	l1 = l + 1
925	l2 = l1 + 1
926	g = d(l)
927	p = (d(l1) - g) / (2.0e0 * e(l))
928	r = pythag(p,1.0e0_rk)
929	d(l) = e(l) / (p + sign(r,p))
930	d(l1) = e(l) * (p + sign(r,p))
931	dl1 = d(l1)
932	h = g - d(l)
933	if (l2 .gt. n) go to 145
934
935	do 140 i = l2, n
936	140 d(i) = d(i) - h
937
938	145 f = f + h
939	! .......... ql transformation ..........
940	p = d(m)
941	c = 1.0e0
942	c2 = c
943	el1 = e(l1)
944	s = 0.0e0
945	mml = m - l
946	! .......... for i=m-1 step -1 until l do -- ..........
947	do 200 ii = 1, mml
948	c3 = c2
949	c2 = c
950	s2 = s
951	i = m - ii
952	g = c * e(i)
953	h = c * p
954	r = pythag(p,e(i))
955	e(i+1) = s * r
956	s = e(i) / r
957	c = p / r
958	p = c * d(i) - s * g
959	d(i+1) = h + s * (c * g + s * d(i))
960	200 continue
961
962	p = -s * s2 * c3 * el1 * e(l) / dl1
963	e(l) = s * p
964	d(l) = c * p
965	tst2 = tst1 + abs(e(l))
966	if (tst2 .gt. tst1) go to 130
967	210 p = d(l) + f
968	! .......... order eigenvalues ..........
969	if (l .eq. 1) go to 250
970	! .......... for i=l step -1 until 2 do -- ..........
971	do 230 ii = 2, l
972	i = l + 2 - ii
973	if (p .ge. d(i-1)) go to 270
974	d(i) = d(i-1)
975	230 continue
976
977	250 i = 1
978	270 d(i) = p
979	290 continue
980
981	go to 1001
982	! .......... set error -- no convergence to an
983	! eigenvalue after 30 iterations ..........
984	1000 ierr = l
985	1001 return
986	end subroutine tql1
987
988	function pythag(a,b) result(pyt)
989	! finds dsqrt(a2+b2) without overflow or destructive underflow
990	implicit none
991	real(kind=rk) :: a,b,pyt
992	real(kind=rk) :: p,r,s,t,u
993
994	p = dmax1(dabs(a),dabs(b))
995	if (p .eq. 0.0d0) go to 20
996	r = (dmin1(dabs(a),dabs(b))/p)**2
997	10 continue
998	t = 4.0d0 + r
999	if (t .eq. 4.0d0) go to 20
1000	s = r/t
1001	u = 1.0d0 + 2.0d0*s
1002	p = u*p
1003	r = (s/u)*2 r
1004	go to 10
1005	20 pyt = p
1006	return
1007	end function pythag
1008
1009
1010	end module pcg_mod
1011
1012

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

Download in other formats:

Original Format