!***********************************************************************
!                              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.                           
!                                                                       
!***********************************************************************
!                                                                       
!***********************************************************************
                                                                        
!-----------------------------------------------------------------------
!  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 
!-----------------------------------------------------------------------
! 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  itercount 
      INTEGER, PARAMETER :: dp = KIND(1.0D0)
                                                                        
      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 :: IA(1000001)!(i*j+1)  row start/end element pointer
      INTEGER :: JA(3990000)!(3*i*j)  column position (absolute)
      REAL(KIND=DP) :: A(3990000)  !( 3*i*j )
      REAL(KIND=DP) :: 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 dpar(length),TMP(1000000,4) 
!-----------------------------------------------------------------------
! Some additional variables to use with the RCI (P)CG solver            
!-----------------------------------------------------------------------
      CHARACTER MATDES(3) 
      DOUBLE PRECISION  solution(1000000) 
      DOUBLE PRECISION expected_sol(1000000) 
      DOUBLE PRECISION  ONE 
      DATA   ONE/1.D0/ 
      DOUBLE PRECISION Euclidean_norm, temp(1000000) 
!-----------------------------------------------------------------------
! 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 initial guess                                          
!-----------------------------------------------------------------------
       DO I=1, N 
         solution(I)=0.D0 
       ENDDO 
       MATDES(1)='D' 
       MATDES(2)='L' 
       MATDES(3)='N' 
      
      WRITE(*, *) ' The system is successfully solved ' 
      WRITE(*, *) ' The following solution obtained ' 
      WRITE(*, *) ' expected_sol solution ' 
  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)

      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 
 1000 CONTINUE 
      STOP 
      END                                           