!=============================================================================== ! Copyright 2004-2018 Intel Corporation All Rights Reserved. ! ! The source code, information and material ("Material") contained herein is ! owned by Intel Corporation or its suppliers or licensors, and title to such ! Material remains with Intel Corporation or its suppliers or licensors. The ! Material contains proprietary information of Intel or its suppliers and ! licensors. The Material is protected by worldwide copyright 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 or other intellectual property rights in the Material ! is granted to or conferred upon you, either expressly, by implication, ! inducement, estoppel or otherwise. Any license under such intellectual ! property rights must be express and approved by Intel in writing. ! ! Unless otherwise agreed by Intel in writing, you may not remove or alter this ! notice or any other notice embedded in Materials by Intel or Intel's ! suppliers or licensors in any way. !=============================================================================== * Content : Intel(R) MKL Cluster Sparse Solver Fortran example * real, double precision, symmetric matrix * ******************************************************************************** C---------------------------------------------------------------------- C Example program to show the use of the "CLUSTER_SPARSE_SOLVER" routine C for solving symmetric linear systems C---------------------------------------------------------------------- program cluster_sparse_solver_sym use mpi implicit none !include 'mkl_cluster_sparse_solver.fi' !include 'mpif.h' C.. Internal solver memory pointer for 64-bit architectures integer*8 pt(64) !TYPE(MKL_CLUSTER_SPARSE_SOLVER_HANDLE) pt(64) C.. All other variables INTEGER*8 maxfct, mnum, mtype, phase, n, nrhs, error, msglvl INTEGER*4 rank, mpi_stat INTEGER*8 iparm(64) INTEGER i integer*8 idum(1) integer*8 nnz integer*8 ii INTEGER j INTEGER*8 , allocatable:: ia(:) INTEGER*8 , allocatable:: ja(:) REAL*8, allocatable:: a(:) REAL*8, allocatable:: b(:) REAL*8, allocatable:: bs(:) REAL*8, allocatable:: x(:) integer*8 i8 REAL*8 ddum(1) REAL*8 res, res0 external mkl_dcsrsymv C.. Fill all arrays containing matrix data. Necessary only on master process DATA nrhs /1/, maxfct /1/, mnum /1/ C.. C.. Initialize MPI. call MPI_INIT(mpi_stat) call MPI_COMM_RANK(MPI_COMM_WORLD, rank, mpi_stat) n=67500 !test fails for this value ! n=65000 !test passes for this value nnz=n*(n-1)/2+n write(*,*)'n=',n write(*,*)'nnz=',nnz !stop allocate(ia(n+1)) allocate(ja(nnz)) allocate(a(nnz)) !stop a=.5d0 ii=0 ia(1)=1 do i=1,n ia(i+1)=ia(i)+n-i+1 do j=i,n ii=ii+1 ja(ii)=j end do a(ii)=1.d0 end do !stop allocate(b(n)) allocate(bs(n)) allocate(x(n)) C.. C.. Set up Cluster Sparse Solver control parameter C.. do i = 1, 64 iparm(i) = 0 enddo iparm(1) = 1 ! no solver default iparm(2) = 2 ! fill-in reordering from METIS iparm(2) = 3 ! OpenMP version of the nested dissection algorithm iparm(6) = 0 ! =0 solution on the first n compoments of x !iparm(8) = 2 ! numbers of iterative refinement steps !iparm(10) = 13 ! perturbe the pivot elements with 1E-13 iparm(11) = 0 ! use nonsymmetric permutation and scaling MPS !iparm(13) = 1 ! maximum weighted matching algorithm is switched-off iparm(40) = 0 ! Input: matrix/rhs/solution stored on master error = 0 ! initialize error flag msglvl = 1 ! print statistical information mtype = -2 ! symmetric, indefinite C.. Initiliaze the internal solver memory pointer. This is only C necessary for the FIRST call of the Cluster Sparse Solver. do i = 1, 64 pt(i)= 0 enddo C.. Reordering and Symbolic Factorization, This step also allocates C all memory that is necessary for the factorization phase = 11 ! only reordering and symbolic factorization call cluster_sparse_solver_64 (pt, maxfct, mnum, mtype, phase, n 1 , a,ia, ja, 1 idum, nrhs, iparm, msglvl, ddum, ddum, MPI_COMM_WORLD, error) if (error .ne. 0) then if (rank.eq.0) write(*,*) 1 'ERROR during symbolic factorization: ', error goto 999 endif if (rank.eq.0) write(*,*) 'Reordering completed ... ' C.. Factorization. phase = 22 ! only factorization call cluster_sparse_solver_64 (pt, maxfct, mnum, mtype, phase, n 1 , a,ia, ja, 1 idum, nrhs, iparm, msglvl, ddum, ddum, MPI_COMM_WORLD, error) if (error .ne. 0) then if (rank.eq.0) write(*,*) ! 'ERROR during numerical factorization: ', error goto 999 endif if (rank.eq.0) write(*,*) 'Factorization completed ... ' C.. Back substitution and iterative refinement phase = 33 ! only solution if (rank.eq.0) then do i = 1, n b(i) = 1.d0 end do endif call cluster_sparse_solver_64 (pt, maxfct, mnum, mtype, phase, n ! , a, ia, ja, 1 idum, nrhs, iparm, msglvl, b, x, MPI_COMM_WORLD, error) if (error.ne.0) then if (rank.eq.0) write(*,*) 'ERROR during during solution: ', ! error goto 999 endif if (rank.eq.0) then write(*,*) 'Solve completed ... ' write(*,*) 'The solution of the system is ' ! do i = 1, n ! write(*,'(" x( ", I1, " ) = ", F19.16)') i, x(i) ! enddo !call mkl_dcsrsymv ('U', n, a, ia, ja, x, bs) bs=0.d0 do i=1,n do i8=ia(i),ia(i+1)-1 j=ja(i8) bs(i)=bs(i)+a(i8)*x(j) if(i.ne.j) bs(j)=bs(j)+a(i8)*x(i) end do end do res = (0.d0,0.d0) res0 = (0.d0,0.d0) do i=1,n res = res + (bs(i)-b(i))**2 res0 = res0 + b(i)**2 enddo print *, 'Relative residual = ', sqrt(abs(res))/ 1 sqrt(abs(res0)) endif C.. Termination and release of memory phase = -1 ! release internal memory call cluster_sparse_solver_64 (pt, maxfct, mnum, mtype, phase 1 , n, ddum, idum, idum, 1 idum, nrhs, iparm, msglvl, ddum, ddum, MPI_COMM_WORLD, error) if (error.ne.0) then if (rank.eq.0) write(*,*) 1 'The following ERROR was detected: ', error goto 999 endif if (rank.eq.0) then if (sqrt(abs(res))/sqrt(abs(res0)).gt.1.d-10) then write(*,*) 'Error: residual is too high!' error = 1 endif if(error.ne.0) then write(*,*) char(10), 'TEST FAILED' else write(*,*) char(10), 'TEST PASSED' endif endif 999 continue call MPI_FINALIZE(mpi_stat) if (error.ne.0.and.rank.eq.0) then stop 1 endif end