!*********************************************************************** ! INTEL CONFIDENTIAL ! Copyright(C) 2005-2008 Intel Corporation. All Rights Reserved. ! The source code contained or described herein and all documents re ! the source code ("Material") are owned by Intel Corporation or its s ! or licensors. Title to the Material remains with Intel Corporatio ! suppliers and licensors. The Material contains trade secrets and pro ! and confidential information of Intel or its suppliers and licens ! Material is protected by worldwide copyright and trade secret ! treaty provisions. No part of the Material may be used, copied, rep ! modified, published, uploaded, posted, transmitted, distributed or d ! in any way without Intel's prior express written permission. ! No license under any patent, copyright, trade secret or other inte ! property right is granted to or conferred upon you by disclosure or ! of the Materials, either expressly, by implication, inducement, est ! otherwise. Any license under such intellectual property rights ! express and approved by Intel in writing. ! !*********************************************************************** ! Content: Intel MKL RCI (P)CG Fortran-77 example ! !*********************************************************************** !----------------------------------------------------------------------- ! Example program for solving symmetric positive definite system of equ ! Full case: full functionality of RCI (P)CG is used. !----------------------------------------------------------------------- PROGRAM rci_pcg use omp_lib IMPLICIT NONE INCLUDE 'mkl_rci.fi' !----------------------------------------------------------------------- ! Define arrays for the upper triangle of the coefficient matrix and ! preconditioner as well as an array for rhs vector ! Compressed sparse row storage is used for sparse representation !----------------------------------------------------------------------- INTEGER RCI_request, itercount INTEGER, PARAMETER :: dp = KIND(1.0D0) ! PARAMETER (N=8) ! DOUBLE PRECISION rhs(N) ! INTEGER IA(9) ! INTEGER JA(18) ! DOUBLE PRECISION A(18) ! Fill all arrays containing matrix data. ! DATA IA /1,5,8,10,12,15,17,18,19/ ! DATA JA & ! & /1, 3, 6,7, & ! & 2,3, 5, & ! & 3, 8, & ! & 4, 7, & ! & 5,6,7, & ! & 6, 8, & ! & 7, & ! & 8/ ! DATA A & ! & /7.D0, 1.D0, 2.D0, 7.D0, & ! & -4.D0, 8.D0, 2.D0, & ! & 1.D0, 5.D0, & ! & 7.D0, 9.D0, & ! & 5.D0, 1.D0, 5.D0, & ! & -1.D0, 5.D0, & ! & 11.D0, & ! & 5.D0/ INTEGER :: ii=100,jj=100,kk=100,i,j,k,ip_exists=1,im_exists=1,jp_exists=1,jm_exists=1,kp_exists=1,n INTEGER nnz,nn,nnzc,nnzc0,nnn,nnzr INTEGER, save :: IA(1000001)!(i*j+1) row start/end element pointer INTEGER, save :: JA(3990000)!(3*i*j) column position (absolute) REAL(KIND=DP), save :: A(3990000) !( 3*i*j ) REAL(KIND=DP), save :: rhs(1000000)!( i*j ) REAL(KIND=DP) :: timest,timeend,time !----------------------------------------------------------------------- ! Allocate storage for the solver ?par and the initial solution vector !----------------------------------------------------------------------- INTEGER length PARAMETER (length=128) INTEGER ipar(length) DOUBLE PRECISION, save :: dpar(length),TMP(1000000,4) !----------------------------------------------------------------------- ! Some additional variables to use with the RCI (P)CG solver !----------------------------------------------------------------------- CHARACTER MATDES(3) DOUBLE PRECISION, save :: solution(1000000) DOUBLE PRECISION, save :: expected_sol(1000000) DOUBLE PRECISION ONE DATA ONE/1.D0/ DOUBLE PRECISION, save :: DNRM2, Euclidean_norm, temp(1000000) EXTERNAL DNRM2 !----------------------------------------------------------------------- ! Initialize the coeffient matrix and expected solution !----------------------------------------------------------------------- n = ii*jj*kk nnz = ii*jj*(kk*3-1) ! expected_sol=1.D0 nnzc = 1 !! nnz accumulative count IA(1)=1 do k = 1,kk do j = 1,jj do i = 1,ii nnzr=0 !!count of nnz(real:include no-path) in a row nn=(k-1)*ii*jj+(j-1)*ii+i JA(nnzc) = nn !!self position nnzc0 = nnzc nnzc = nnzc+1 JA(nnzc) = nn+1 !!normal x_p position if (i/=ii) then if(ip_exists==1)then A(nnzc) = -1d0 !!normal x_p value nnzr=nnzr+1 else A(nnzc) = 0d0 endif nnzc = nnzc+1 endif if (i==1) then JA(nnzc) = nn+ii-1 if(im_exists==1)then !! x-periodic material exist A(nnzc) =-1d0 ! -1d0 !!... value nnzr=nnzr+1 else A(nnzc) = 0d0 endif nnzc = nnzc+1 endif if (j/=jj) then JA(nnzc) = nn+ii !!normal y_p position if(jp_exists==1)then A(nnzc) = -1d0 !!... value nnzr=nnzr+1 else A(nnzc) = 0d0 endif nnzc = nnzc+1 endif if (j==1) then JA(nnzc) = nn+(jj-1)*ii !!periodic y_m position if(jm_exists==1)then A(nnzc) = -1d0 !!... value nnzr=nnzr+1 else A(nnzc) = 0d0 endif nnzc = nnzc+1 endif if (k/=kk) then JA(nnzc) = nn+jj*ii !!normal z_p position if(kp_exists==1)then A(nnzc) = -1d0 !!... value nnzr=nnzr+1 else A(nnzc) = 0d0 endif nnzc = nnzc+1 A(nnzc0) =6d0 !dble(nnzr) !!self value else A(nnzc0) =5d0 !dble(nnzr) !!self value endif IA(nn+1)=nnzc enddo enddo enddo rhs=0d0 do k = 40,60 do j = 20,40 do i = 10,60 nn=(k-1)*ii*jj+(j-1)*ii+i rhs(nn) = 1.d-5 enddo end do enddo !----------------------------------------------------------------------- ! Initialize the right hand side through matrix-vector product !----------------------------------------------------------------------- ! CALL MKL_DCSRSYMV('U', N, A, IA, JA, expected_sol, rhs) !----------------------------------------------------------------------- ! Initialize the initial guess !----------------------------------------------------------------------- DO I=1, N solution(I)=0.D0 ENDDO MATDES(1)='D' MATDES(2)='L' MATDES(3)='N' !----------------------------------------------------------------------- ! Initialize the solver !----------------------------------------------------------------------- ! timest = omp_get_wtime() CALL dcg_init(N, solution,rhs, RCI_request,ipar,dpar,TMP) IF (RCI_request .NE. 0 ) GOTO 999 !----------------------------------------------------------------------- ! Set the desired parameters: ! INTEGER parameters: ! set the maximal number of iterations to 100 ! LOGICAL parameters: ! run the Preconditioned version of RCI (P)CG with preconditioner C_inve ! DOUBLE PRECISION parameters ! - !----------------------------------------------------------------------- ipar(5)=1000 ipar(11)=1 !----------------------------------------------------------------------- ! Check the correctness and consistency of the newly set parameters !----------------------------------------------------------------------- CALL dcg_check(N,solution,rhs,RCI_request,ipar,dpar,TMP) IF (RCI_request .NE. 0 ) GOTO 999 !----------------------------------------------------------------------- ! Compute the solution by RCI (P)CG solver ! Reverse Communications starts here !----------------------------------------------------------------------- 1 CALL dcg(N,solution,rhs,RCI_request,ipar,dpar,TMP) !----------------------------------------------------------------------- ! If RCI_request=0, then the solution was found according to the request ! stopping tests. In this case, this means that it was found after 100 ! iterations. !----------------------------------------------------------------------- ! print*,'RCI=',RCI_request IF (RCI_request .EQ. 0) THEN GOTO 700 !----------------------------------------------------------------------- ! If RCI_request=1, then compute the vector A*TMP(:,1) ! and put the result in vector TMP(:,2) !----------------------------------------------------------------------- ELSEIF (RCI_request .EQ. 1) THEN CALL MKL_DCSRSYMV('U', N, A, IA, JA, TMP, TMP(1,2)) GOTO 1 !----------------------------------------------------------------------- ! If RCI_request=2, then do the user-defined stopping test: compute the ! Euclidean norm of the actual residual using MKL routines and check if ! it is less than 1.D-8 !----------------------------------------------------------------------- ELSEIF (RCI_request .EQ. 2) THEN CALL MKL_DCSRSYMV('U', N, A, IA, JA, solution, temp) CALL DAXPY(N,-1.D0,rhs,1,temp,1) Euclidean_norm = DNRM2(N,temp,1) IF (Euclidean_norm .GT. 1.D-8) THEN !----------------------------------------------------------------------- ! The solution has not been found yet according to the user-defined stop ! test. Continue RCI (P)CG iterations. !----------------------------------------------------------------------- GOTO 1 ELSE !----------------------------------------------------------------------- ! The solution has been found according to the user-defined stopping tes !----------------------------------------------------------------------- GOTO 700 END IF !----------------------------------------------------------------------- ! If RCI_request=3, then compute apply the preconditioner matrix C_inver ! on vector TMP(:,3) and put the result in vector TMP(:,4) !----------------------------------------------------------------------- ELSEIF (RCI_request .EQ. 3) THEN CALL MKL_DCSRSV('N', N, ONE, MATDES, & & A, JA, IA, IA(2), TMP(1,3), TMP(1,4)) GOTO 1 ELSE !----------------------------------------------------------------------- ! If RCI_request=anything else, then dcg subroutine failed ! to compute the solution vector: solution(N) !----------------------------------------------------------------------- GOTO 999 ENDIF !----------------------------------------------------------------------- ! Reverse Communication ends here ! Get the current iteration number !----------------------------------------------------------------------- 700 CALL dcg_get(N,solution,rhs,RCI_request,ipar,dpar,TMP, & & itercount) ! print*,'RCI=',RCI_request !----------------------------------------------------------------------- ! Print solution vector: solution(N) and number of iterations: itercount !----------------------------------------------------------------------- ! timeend = omp_get_wtime() ! time = time + (timest - timeend) ! !$omp parallel do& ! !$omp& private(i) ! do i= 1,36 ! WRITE(*, *) i,time,omp_get_thread_num() ! end do ! !$omp end parallel do WRITE(*, *) ' The system is successfully solved ' WRITE(*, *) ' The following solution obtained ' ! WRITE(*,800) (solution(i),i =1,N) WRITE(*, *) ' expected_sol solution ' ! WRITE(*,800)(expected_sol(i),i =1,N) 800 FORMAT(4(F10.3)) WRITE(*,900)(itercount) 900 FORMAT(' Number of iterations: ',1(I4)) open(100,file='testcg2.dat') DO i = 1, n WRITE(100,*) solution(i) END DO close(100) CALL MKL_DCSRSYMV('U', N, A, IA, JA, solution, expected_sol) open(100,file='Axcg.dat') DO i = 1, n WRITE(100,*) expected_sol(i) END DO close(100) open(100,file='rhscg.dat') DO i = 1, n WRITE(100,*) rhs(i) END DO close(100) GOTO 1000 999 WRITE(*,*) 'Solver returned error code ', RCI_request STOP 1000 CONTINUE STOP END