Context Navigation

pcg.F90 @ b378423

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

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

Download in other formats:

Original Format