Context Navigation

pcg.F90 @ 8e582be

Revision 8e582be, 33.9 KB checked in by chris <chris@…>, 7 years ago (diff)

Changed license header

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

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

Download in other formats:

Original Format