From 2cabc9f82bddb0ce8e4e226c8ffa6a0d0446bdd6 Mon Sep 17 00:00:00 2001 From: Jan Mandel Date: Tue, 7 Sep 2021 19:28:19 -0600 Subject: [PATCH] adding prints, computing residual in matri-vector multiply #4 --- femwind/fortran/module_multigrid.f90 | 10 ++++++---- femwind/fortran/module_ndt_mult.f90 | 13 ++++++++++--- femwind/fortran/module_sweeps.f90 | 6 +++--- femwind/fortran/ndt_mult_test.f90 | 3 ++- femwind/fortran/sweeps_test.f90 | 4 ++-- 5 files changed, 23 insertions(+), 13 deletions(-) diff --git a/femwind/fortran/module_multigrid.f90 b/femwind/fortran/module_multigrid.f90 index 409e6c1..7fc9d07 100644 --- a/femwind/fortran/module_multigrid.f90 +++ b/femwind/fortran/module_multigrid.f90 @@ -181,7 +181,7 @@ integer:: & integer::it, maxit, nit, nsmooth, cycles logical::coarse,coarsest -real:: restol, norm_f, norm_res, diftol, reldif +real:: restol, norm_f, norm_res, diftol, siz, reldif, relres, rhsiz character(len=10)::it_kind !*** executable @@ -254,7 +254,8 @@ do it=1,maxit ifms, ifme, kfms, kfme, jfms, jfme, & ! fire memory bounds ifps, ifpe, kfps, kfpe, jfps, jfpe, & ! fire patch bounds ifts, ifte, kfts, kfte, jfts,jfte, & ! fire tile bounds - mg(l)%Kglo, mg(l)%lambda, mg(l)%f, mg(l)%res) + mg(l)%Kglo, mg(l)%lambda, mg(l)%f, mg(l)%res, rhsiz, relres) + if(params%print_level>=l)print *,'level=',l,' it=',it,'mult rhsiz=',siz,' rel res =',reldif if(params%debug_level >=l)call write_array(mg(l)%res(ifts: ifte, kfts: kfte, jfts:jfte),'cc_res') ! restriction: f_coarse = R*residual call restriction( & @@ -300,7 +301,8 @@ do it=1,maxit ifms, ifme, kfms, kfme, jfms, jfme, & ! memory dimensions ifps, ifpe, kfps, kfpe, jfps, jfpe, & ! fire patch bounds ifts, ifte, kfts, kfte, jfts, jfte, & ! tile dimensions - mg(l)%Kglo, mg(l)%f, mg(l)%lambda, reldif) + mg(l)%Kglo, mg(l)%f, mg(l)%lambda, siz, reldif) + if(params%print_level>=l)print *,'level=',l,' it=',it,'sweeps siz=',siz,' rel diff=',reldif if(reldif < diftol)then if(params%print_level>=l)print *,'level ',l,' iteration ',it,' rel diff=',reldif,'< diftol=',diftol,' exiting' exit @@ -315,7 +317,7 @@ do it=1,maxit ifms, ifme, kfms, kfme, jfms, jfme, & ! fire memory bounds ifps, ifpe, kfps, kfpe, jfps, jfpe, & ! fire patch bounds ifts, ifte, kfts, kfte, jfts,jfte, & ! fire tile bounds - mg(l)%Kglo, mg(l)%lambda, mg(l)%f, mg(l)%res) + mg(l)%Kglo, mg(l)%lambda, mg(l)%f, mg(l)%res, siz, relres) if(params%debug_level >=l)call write_array(mg(l+1)%res(ifcts: ifcte, kfcts: kfcte, jfcts:jfcte),'it_res') norm_res = norm2( mg(l)%res(ifts: ifte, kfts: kfte, jfts: jfte)) ! residual norm if (mod(it,params%nsmooth+1)==params%nsmooth)then diff --git a/femwind/fortran/module_ndt_mult.f90 b/femwind/fortran/module_ndt_mult.f90 index 66dd4a9..ab4d7aa 100644 --- a/femwind/fortran/module_ndt_mult.f90 +++ b/femwind/fortran/module_ndt_mult.f90 @@ -9,7 +9,7 @@ subroutine ndt_mult( & ifms, ifme, kfms,kfme, jfms, jfme, & ifps, ifpe, kfps, kfpe, jfps, jfpe, & ! fire patch bounds ifts, ifte, kfts, kfte, jfts,jfte, & - kmat, lambda, y, r) + kmat, lambda, y, r, siz, relres) implicit none @@ -27,17 +27,19 @@ integer, parameter:: msize = 14 real, intent(in), dimension(ifms:ifme,kfms:kfme,jfms:jfme,msize):: kmat ! global stiffness matrix real,intent(in), dimension(ifms:ifme,kfms:kfme,jfms:jfme):: lambda, y ! input vectors real,intent(out), dimension(ifms:ifme,kfms:kfme,jfms:jfme):: r ! output vector +real, intent(out):: siz, relres !*** local integer:: i,j,k,ie,je,ke -real:: t +real:: t, res !** executable ie = snode(ifte,ifde,+1) je = snode(jfte,jfde,+1) ke = snode(kfte,kfde,+1) - +siz = 0. +res = 0. do j=jfts,je do k=kfts,ke @@ -71,9 +73,14 @@ do j=jfts,je kmat(i ,k ,j ,13)*lambda(i ,k+1,j+1) + & kmat(i ,k ,j ,14)*lambda(i+1,k+1,j+1) ) r(i,k,j) = t + siz = max(siz,abs(y(i,k,j))) + res = max(res,abs(t)) enddo enddo enddo + +relres = res/max(tiny(siz),siz) + end subroutine ndt_mult diff --git a/femwind/fortran/module_sweeps.f90 b/femwind/fortran/module_sweeps.f90 index 27186a6..f5c8b86 100644 --- a/femwind/fortran/module_sweeps.f90 +++ b/femwind/fortran/module_sweeps.f90 @@ -7,7 +7,7 @@ subroutine sweeps(ifds, ifde, kfds,kfde, jfds, jfde, & ! fire ifms, ifme, kfms,kfme, jfms, jfme, & ifps, ifpe, kfps, kfpe, jfps, jfpe, & ! fire patch bounds ifts, ifte, kfts, kfte, jfts,jfte, & - Kmat, F, x, reldif) !input and output matrices/vectors + Kmat, F, x, siz, reldif) !input and output matrices/vectors implicit none @@ -24,12 +24,12 @@ integer, parameter:: msize = 14 real, intent(in), dimension(ifms:ifme,kfms:kfme,jfms:jfme,msize):: Kmat ! global stiffness matrix real, intent(in), dimension(ifms:ifme,kfms:kfme,jfms:jfme):: F ! global load vector real,intent(out), dimension(ifms:ifme,kfms:kfme,jfms:jfme):: x ! output vector -real, intent(out):: reldif ! relative difference betweet input and output, in max norm +real, intent(out):: siz, reldif ! size of input, relative difference betweet input and output, in max norm !*** local integer:: i, j, k, r1, r2, ie,je, ke -real:: t, dif, siz, old, new +real:: t, dif, old, new !*** executable diff --git a/femwind/fortran/ndt_mult_test.f90 b/femwind/fortran/ndt_mult_test.f90 index eca9f3a..e4e8bf8 100644 --- a/femwind/fortran/ndt_mult_test.f90 +++ b/femwind/fortran/ndt_mult_test.f90 @@ -16,6 +16,7 @@ integer :: msize, & ifps, ifpe, kfps, kfpe, jfps, jfpe, & ! fire patch bounds ifts, ifte, kfts, kfte, jfts,jfte ! fire tile bounds integer :: i,j,k,jx +real:: siz,relres ! read input arrays in ikj index ordering and tight bounds call read_array_nd(a,s,'kmat') @@ -77,7 +78,7 @@ call ndt_mult( & ifms, ifme, kfms, kfme, jfms, jfme, & ! fire memory bounds ifps, ifpe, kfps, kfpe, jfps, jfpe, & ! fire patch bounds ifts, ifte, kfts, kfte, jfts,jfte, & ! fire tile bounds - kmat, u, y, r ) + kmat, u, y, r, siz, relres) write(*,'(a,3i8)')'copying -y to array size ',n allocate(y_m(1:n(1),1:n(2),1:n(3))) diff --git a/femwind/fortran/sweeps_test.f90 b/femwind/fortran/sweeps_test.f90 index 88a3cb4..989afb7 100644 --- a/femwind/fortran/sweeps_test.f90 +++ b/femwind/fortran/sweeps_test.f90 @@ -9,7 +9,7 @@ real, pointer:: Kmat(:,:,:,:), x(:,:,:),F(:,:,:), & ! fortran is not case sensi Kmat_m(:,:,:,:), x_m(:,:,:), F_m(:,:,:) real, pointer::a(:) integer :: s(4),n(3) -real::reldif +real::reldif, siz integer :: msize, & ifds, ifde, kfds, kfde, jfds, jfde, & ! fire domain bounds @@ -74,7 +74,7 @@ call sweeps( & ifms, ifme, kfms, kfme, jfms, jfme, & ! fire memory bounds ifps, ifpe, kfps, kfpe, jfps, jfpe, & ! fire patch bounds ifts, ifte, kfts, kfte, jfts,jfte, & ! fire tile bounds - Kmat, F, x, reldif) + Kmat, F, x, siz, reldif) write(*,'(a,3i8)')'copying the output data to array size ',n allocate(x_m(1:n(1),1:n(2),1:n(3))) -- 2.11.4.GIT