Conjugate Gradient source Code Example

[Juan_Diego_E_]

Good afternoon,

I would like to have some help in the process of linking MKL with Visual Studio in Fortran Language. I am trying to run the example "cg_jacobi_precon" that comes in the Intel Folder, but I have not succeed. I am actually using F90. Can I have any advice, please, on which source codes do I need to include in my project as well as header files, and if I have to set any address in the Fortran and/or Linker properties of the project?

Please find below the actual example (I guess it is in F77) from Intel.

Thank you.

Regards,

JD.

!*******************************************************************************
!                              INTEL CONFIDENTIAL
!   Copyright(C) 2005-2010 Intel Corporation. All Rights Reserved.
!   The source code contained  or  described herein and all documents related to
!   the source code ("Material") are owned by Intel Corporation or its suppliers
!   or licensors.  Title to the  Material remains with  Intel Corporation or its
!   suppliers and licensors. The Material contains trade secrets and proprietary
!   and  confidential  information of  Intel or its suppliers and licensors. The
!   Material  is  protected  by  worldwide  copyright  and trade secret laws and
!   treaty  provisions. No part of the Material may be used, copied, reproduced,
!   modified, published, uploaded, posted, transmitted, distributed or disclosed
!   in any way without Intel's prior express written permission.
!   No license  under any  patent, copyright, trade secret or other intellectual
!   property right is granted to or conferred upon you by disclosure or delivery
!   of the Materials,  either expressly, by implication, inducement, estoppel or
!   otherwise.  Any  license  under  such  intellectual property  rights must be
!   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 equations.
!  Full case: full functionality of RCI (P)CG is used.
!---------------------------------------------------------------------------
      PROGRAM rci_pcg_f77_test_3
      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 N, RCI_request, itercount, expected_itercount, i
      PARAMETER (N=8)
      PARAMETER (expected_itercount=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 /1,  3,     6,7,
     2    2,3,   5,
     3      3,         8,
     4         4,    7,
     5           5,6,7,
     6             6,  8,
     7               7,
     8                 8/
      DATA A
     1 /7.D0,      1.D0,           2.D0, 7.D0,
     2       -4.D0, 8.D0,     2.D0,
     3             1.D0,                      5.D0,
     4                  7.D0,      9.D0,
     5                       5.D0, 1.D0, 5.D0,
     6                            -1.D0,      5.D0,
     7                                  11.D0,
     8                                        5.D0/

!---------------------------------------------------------------------------
! Allocate storage for the solver ?par and the initial solution vector
!---------------------------------------------------------------------------
      INTEGER length
      PARAMETER (length=128)
      INTEGER ipar(length)
      DOUBLE PRECISION dpar(length),TMP(N,4)
!---------------------------------------------------------------------------
! Some additional variables to use with the RCI (P)CG solver
!---------------------------------------------------------------------------
      CHARACTER MATDES(3)
      DOUBLE PRECISION  solution(N)
      DOUBLE PRECISION expected_sol(N)
      DATA expected_sol/1.D0, 0.D0, 1.D0, 0.D0, 1.D0, 0.D0, 1.D0, 0.D0/
      DOUBLE PRECISION  ONE
      DATA   ONE/1.D0/
      DOUBLE PRECISION DNRM2, Euclidean_norm, temp(N)
      EXTERNAL DNRM2
!---------------------------------------------------------------------------
! 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
!---------------------------------------------------------------------------
      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_inverse
! DOUBLE PRECISION parameters
! -
!---------------------------------------------------------------------------
      ipar(5)=100
      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 requested
! stopping tests. In this case, this means that it was found after 100
! iterations.
!---------------------------------------------------------------------------
      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 stopping
! test. Continue RCI (P)CG iterations.
!---------------------------------------------------------------------------
          GOTO 1
        ELSE
!---------------------------------------------------------------------------
! The solution has been found according to the user-defined stopping test
!---------------------------------------------------------------------------
          GOTO 700
        END IF
!---------------------------------------------------------------------------
! If RCI_request=3, then compute apply the preconditioner matrix C_inverse
! 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 solution vector: solution(N) and number of iterations: itercount
!---------------------------------------------------------------------------
      WRITE(*, *) ' The system has been solved '
      WRITE(*, *) ' The following solution obtained '
      WRITE(*,800) (solution(i),i =1,N)
      WRITE(*, *) ' expected solution '
      WRITE(*,800)(expected_sol(i),i =1,N)
800   FORMAT(4(F10.3))
      WRITE(*,900)(itercount)
900   FORMAT(' Number of iterations: ',1(I2))
      DO I=1,N
         expected_sol(I)=expected_sol(I)-solution(I)
      ENDDO

      Euclidean_norm=DNRM2(N,expected_sol,1)

!---------------------------------------------------------------------------
! Release internal MKL memory that might be used for computations
! NOTE: It is important to call the routine below to avoid memory leaks
! unless you disable MKL Memory Manager
!---------------------------------------------------------------------------
      CALL MKL_FREEBUFFERS

      IF (itercount.EQ.expected_itercount .AND.
     1                                   Euclidean_norm.LE.1.0D-12) THEN
         WRITE( *,'(A,A)') 'This example has successfully PASSED',
     1 ' through all steps of computation!'
         STOP 0
      ELSE
         WRITE( *,'(A,A,A,I5,A,A,A,E12.5,A)') 'This example may have',
     1 ' FAILED as either the number of iterations differs from the',
     2 ' expected number of iterations ',expected_itercount,' or the',
     3 ' computed solution differs much from the expected solution (',
     4 'Euclidean norm is ',Euclidean_norm,'), or both.'
         STOP 1
      ENDIF
!---------------------------------------------------------------------------
! Release internal MKL memory that might be used for computations
! NOTE: It is important to call the routine below to avoid memory leaks
! unless you disable MKL Memory Manager
!---------------------------------------------------------------------------
999   WRITE( *,'(A,A)') 'This example FAILED as the solver has',
     1 ' returned the ERROR code', RCI_request
      CALL MKL_FREEBUFFERS
      STOP 1

      END

[Juan_Diego_E_]

Good afternoon,

I just wanted to clarify that I managed to run the code (I had a typo in one of the lines).

I would also like to ask out there, if FGMRES MKL subroutine is threaded or not? (let say the example that comes with the ILU0 preconditioner, how can I run it in parallel?). If it is not threaded, I would like to have any advice on how can I use OpenMP or any compiler options (from Visual Studio Properties) to parallelize this FGMRES subroutine. I am solving a finite difference problem with multiple time steps, therefore I have to call a lot of times a subroutine that itself calls DFMGRES (like in the script above).

Thank you very much in advance.

Regards,

J.D.

[Gennady_F_Intel]

both Preconditioners  (ILUT and IL0) from MKL are not threaded