program CMRH !*************************************************************** ! CMRH program !*************************************************************** external SSWAP,DCOPY,DSCAL,DGEMV,DSWAP intrinsic dabs,dsqrt !******************************************** ! VARIABLES !******************************************** integer ::n,code,i,j,k,choix,solchoice,iprint real(kind=8) ::init(5) real(kind=8),allocatable,dimension(:,:) ::A real(kind=8) ::btemp,bmax,btemp2,err,res real(kind=8),allocatable,dimension(:) ::b,sol,rs,temp,x,work,tempn,c,s,temp2,xx,yy,valeurs integer,allocatable,dimension(:) ::p real*8 ::dnorm2,D_inv,beta,dnorm,dnrm2,tol,alpha,H_k1,dgam,r integer ::indice,idamax,k1,n_k real dtime, delta, delta0, da(2), dt1, dt2, dt3 real etime, eelta, eelta0, ea(2), et1, et2, et3 real t1, t2, t3, tm1, tm2, tm3 real dtime1, delta1, delta01, da1(2) real etime1, eelta1, eelta01, ea1(2),l,p1,p2 integer nb_ligne real ::rand da(1)=0.d0 da(2)=1.d0 !call itime(timearray) ! Get the current time !i = rand ( timearray(1)+timearray(2)+timearray(3) ) !************************************ ! Matrix Initialisation !************************************ open(unit=10,file='inputfile.dat',status='old') read(10,*)init close(10) iprint=int(init(1)) n=int(init(2)) choix=int(init(3)) solchoice=int(init(4)) tol=init(5) !call init(iprint,n,choix,solchoice) allocate(b(n),rs(n),x(n),A(n,n),work(n),sol(n)) if (choix<3)then call matrice(A,n,n,2.d0,1.d0,choix) else allocate(valeurs(n)) open(unit=30,file='matrixfile.dat',status='old') do i=1,n read(30,*)valeurs A(i,:)=valeurs(:) enddo close(30) endif !******************************************** ! End Matrix initialisation !******************************************** !******************************************** ! Solution initialisation !******************************************** call solution(A,b,sol,n,n,solchoice) !******************************************** ! End solution initialisation !******************************************** !******************************************** ! Permutation Vector initialisation ! and computation of initial residual !******************************************** allocate(p(n)) do i=1,n p(i)=i end do dnorm2=dnrm2(n,b,1) if (iprint==2)then print*,'' print*,'||r0|| =',dnorm2 print*,'' end if !******************************************** ! Computation of max b(i) !******************************************** indice = idamax(n,b,1) beta = b(indice) D_inv = 1.d0/beta !******************************************** ! First permutation !******************************************** if (indice .ne.1) then call dswap(1,b(indice),1,b(1),1) call dswap(n,A(1,indice),1,A(1,1),1) call dswap(n,A(indice,1),n,A(1,1),n) call sswap(1,p(indice),1,p(1),1) end if !******************************************** ! End permutation !******************************************** ! call dscal(n,0.d0,x,1) call dscal(n,D_inv,b,1) rs(1)=beta dnorm=dabs(beta) allocate(c(n),s(n)) call dscal(n,0.d0,x,1) ! Fin initialisation de A !print*,'Fin initialisation des matrices' k=0 if (iprint==2)then write(*,*)'iter., pseudo-residual norm' end if !tol=1.d-7 delta0 = dtime(da) eelta0 = etime(ea) open(unit=30,file='residual.dat',status='unknown') write(30,*)'Iteration, Pseudo-Residual Norm' !******************************************** ! Loop !******************************************** do while((dnorm.gt.tol).and.(k.lt.n)) k=k+1; n_k=n-k; k1=k+1; if (iprint==2)then write(*,10)k-1,dnorm end if write(30,10)k-1,dnorm 10 format(i5,2x,e12.4) !******************************************** ! workl = A*b !******************************************** call dcopy(n,A(1,k),1,work,1) call dgemv('N',n,n_k,1.d0,A(1,k1),n,b(k1),1,1.d0,work,1) ! ! call dscal(nl,0.d0,workl,1) ! ! do j=k1,n ! ! do i=1,nl ! ! workl(i)=workl(i)+AL(i,j)*b(j) ! ! end do ! ! end do !******************************************** ! A(:,k)=b !******************************************** call dcopy(n_k,b(k1),1,A(k1,k),1) !******************************************** ! j Loop !******************************************** do j=1,k !******************************************** ! A(j,k)=workl(j) ! workl(j)=0 !******************************************** A(j,k)=work(j) work(j)=0.d0 alpha=-A(j,k) !******************************************** ! workl=workl - alpha A(:,j) !******************************************** call daxpy(n-j,alpha,A(j+1,j),1,work(j+1),1) end do !******************************************** ! End j Loop !******************************************** if (k.lt.n)then !******************************************** ! Search indice i_0 !******************************************** indice=k+idamax(n_k,work(k1),1) H_k1=work(indice) !******************************************** ! Permutation !******************************************** if (indice.ne.k1) then call sswap(1,p(indice),1,p(k1),1) call dswap(1,work(indice),1,work(k1),1) call dswap(n,A(1,indice),1,A(1,k1),1) call dswap(n,A(indice,1),n,A(k1,1),n) end if !******************************************** ! End Permutation !******************************************** if (H_k1.ne.0.d0)then !******************************************** ! b=b/(b)_(i_0) !******************************************** D_inv=1.d0/H_k1 call dcopy(n,work,1,b,1) call dscal(n,D_inv,b,1) end if else H_k1=0.d0 call dscal(n,0.d0,b,1) end if !******************************************** ! Application of Givens rotations !******************************************** if (k.gt.1)then do i=2,k btemp=A(i-1,k) A(i-1,k)=c(i-1)*btemp+s(i-1)*A(i,k) A(i,k)=-s(i-1)*btemp+c(i-1)*A(i,k) end do end if if (k.le.n)then btemp=A(k,k) dgam=dsqrt(btemp**2+H_k1**2) c(k)=dabs(btemp)/dgam s(k)=(H_k1*dabs(btemp))/(btemp*dgam) rs(k1)=-s(k)*rs(k) rs(k)=c(k)*rs(k) dnorm=dabs(rs(k1)) A(k,k)=c(k)*btemp+s(k)*H_k1 end if !******************************************** ! End Givens rotations !******************************************** end do !******************************************** ! Solve H_k d_k = alpha e_1 !******************************************** call dtrsv('U','N','N',k,A,n,rs,1) call dcopy(k,rs,1,work,1) call dtrmv('L','N','U',k,A,n,rs,1) call dgemv('N',n-k,k,1.d0,A(k1,1),n,work,1,0.d0,c,1) call daxpy(k,1.d0,rs,1,x,1) call daxpy(n-k,1.d0,c,1,x(k1),1) do j=1,n b(p(j))=x(j) end do !!*********************************************** !!allocate(b(n),rs(n),x(n),A(n,n),work(n),sol(n)) !if (choix<3)then !call matrice(A,n,n,2.d0,1.d0,choix) !else !allocate(valeurs(n)) !open(unit=30,file='matrixfile.dat',status='old') !do i=1,n !read(30,*)valeurs !A(i,:)=valeurs(:) !enddo !close(30) !endif !!******************************************** !! End Matrix initialisation !!******************************************** !!******************************************** !! Solution initialisation !!******************************************** !call solution(A,work,sol,n,n,solchoice) !!********************************************* !! Residual and Error Norm !!********************************************* !call dgemv('N',n,n,-1.d0,A(k1,1),n,b,1,1.d0,work,1) !dnorm=dnorm2(n,work,1) !call dcopy(n,sol,1,rs,1) !call daxpy(n,-1.d0,b,1,rs,1) !dnorm2=dnorm2(n,rs,1) !!*********************************************** delta = dtime(da) eelta = etime(ea) dt2 = delta-delta0 et2 = eelta-eelta0 t2 = dt2+et2 tm2 = 0.5*(t2) !tm2=0.d0 if (iprint >= 2) then print*,'Results : iter, pseudo-residual norm, time' write(*,40)k,dnorm, tm2 end if write(30,20)k,dnorm close(30) open(unit=100,file='tempsseq.dat',status='unknown',access='append') write(unit=100,FMT=*)'Iteration Number :',k,' temps :',tm2 close(100) !20 format(i5,2x,e12.4,e12.4) if (iprint==2) then write(*,*) write(*,*)' Swap, CMRH Solution, Exact Solution : ' write(*,*)' ______________________________________________ ' if (n.lt.10)then do i=1,n write(*,20)p(i), b(i), sol(i) end do else do i=1,10 write(*,20)p(i), b(i), sol(i) end do end if end if 30 format(A37,4(f14.7,1x)) 20 format(1x,i5,2(1x,f14.8)) 40 format(1x,i5,2(1x,e12.4)) deallocate(A,b,sol,work,c,s,rs,x) deallocate(p) end program CMRH !******************************************** ! End CMRH !******************************************** !******************************************** ! Matrix Choice !******************************************** subroutine matrice(A,nl,n,alpha,beta,choix) integer ::n,i,j,choix,nl real(kind=8),dimension(nl,n)::A real*8 ::alpha,beta if (choix==1) then do i = 1, nl do j = 1, i A(i,j) = (2.d0*dble(j)-1.d0)/dble(n-i+j) end do do j = i+1, n A(i,j) = (2.d0*dble(i)-1.d0)/dble(n-i+j) end do end do end if if (choix==2)then do i=1,nl do j=1,i A(i,j)=dble(i) end do do j=i+1,n A(i,j)=dble(j) end do end do end if end subroutine matrice !******************************************** ! Solution Choice !******************************************** subroutine solution(A,b,sol,nl,n,solchoice) integer ::n,i,nl,solchoice real(kind=8),dimension(n,n) ::A real(kind=8),dimension(n) ::sol real(kind=8),dimension(nl) ::b do i=1,n if (solchoice==1) then sol(i)=1.d0 end if if (solchoice==2) then sol(i)=dble(n-i+1) end if if (solchoice==3)then sol(i)=rand(0) end if end do call dgemv('N',nl,n,1.d0,A,nl,sol,1,0.d0,b,1) end subroutine solution !subroutine init(iprint,n,choix,choixsol) !integer :: iprint,n,choix,choixsol !print*,'Would you print results ? : 1-No / 2-All / 3-Time and final residual norm' !read*,iprint !print*,'Enter matrix size : ' !read*,n !print*,'Enter matrix choice : (more details in README.txt) ' !print*,'1- A(i,j)=2j-1 / n-i+j for j=1,k ' !print*,' 2i-1 / n-i+j for j=i+1,n ' !print*,'2- A(i,j)=i for j=1,i ' !print*,' j for j=i+1,n ' !print*,'3- Matrix File (copy in matrixfile.dat) ' !read*,choix !print*,'Enter solution choice : 1- [1,...,1]^T / 2- [n,...,1]^T ' !read*,choixsol !endsubroutine init