!
! This sample program shows how to build and solve a sparse linear
!
! The program  solves a linear system based on the partial differential
! equation 
!
! 
!
!     the equation generated is:
!   b1 d d (u)  b2  d d (u)    a1 d (u))  a2 d (u)))   
! -   ------   -    ------  +    ----- +  ------     + a3 u = 0
!      dx dx         dy dy         dx        dy        
!
! 
! with  Dirichlet boundary conditions on the unit cube 
!
!    0<=x,y<=1
! 
! The equation is discretized with finite differences and uniform stepsize;
! the resulting  discrete  equation is
!
! ( u(x,y)(2b1+2b2+a1+a2)+u(x-1,y)(-b1-a1)+u(x,y-1)(-b2-a2)+
!  -u(x+1,y)b1-u(x,y+1)b2)*(1/h**2)
!
! Example taken from: C.T.Kelley
!    Iterative Methods for Linear and Nonlinear Equations
!    SIAM 1995
!
!
! In this sample program the index space of the discretized
! computational domain is first numbered sequentially in a standard way, 
! then the corresponding vector is distributed according to an HPF BLOCK
! distribution directive.
!
! Boundary conditions are set in a very simple way, by adding 
! equations of the form
!
!   u(x,y) = rhs(x,y)
!
Program PP2D
  USE F90SPARSE
  Implicit none

  interface 
    !.....user passed subroutine.....
    subroutine part_block(glob_index,n,np,pv,nv)
      INTEGER, INTENT(IN)  :: GLOB_INDEX, N, NP
      INTEGER, INTENT(OUT) :: NV
      INTEGER, INTENT(OUT) :: PV(*) 
    end subroutine part_block
  end interface
  ! input parameters
  Character :: CMETHD*10, PREC*10, AFMT*5
  Integer      :: IDIM, IRET

  ! Miscellaneous 
  Integer, Parameter   :: IZERO=0, IONE=1
  Character, PARAMETER :: ORDER='R'
  INTEGER              :: IARGC,CONVERT_DESCR,dim, CHECK_DESCR
  REAL(KIND(1.D0)), PARAMETER :: DZERO = 0.D0, ONE = 1.D0
  REAL(KIND(1.D0)) :: MPI_WTIME, T1, T2, TPREC, TSOLVE, T3, T4 
  EXTERNAL  MPI_WTIME

  ! Sparse Matrix and preconditioner
  TYPE(D_SPMAT) :: A,  L, U, H
  TYPE(D_PREC) :: PRE
  ! Descriptor
  TYPE(desc_type)    :: DESC_A, DESC_A_OUT
  ! Dense Matrices
  REAL(KIND(1.d0)), POINTER :: B(:), X(:), D(:),LD(:)
  INTEGER, pointer :: WORK(:)
  ! BLACS parameters
  INTEGER            :: nprow, npcol, icontxt, iam, np, myprow, mypcol
  
  ! Solver parameters
  INTEGER            :: ITER, ITMAX,IERR,ITRACE, METHD,IPREC, ISTOPC,&
       & IPARM(20), ML
  REAL(KIND(1.D0))   :: ERR, EPS, RPARM(20)
   
  ! Other variables
  INTEGER            :: I,INFO
  INTEGER            :: INTERNAL, M,II
 
  ! Initialize BLACS  
  CALL BLACS_PINFO(IAM, NP)
  CALL BLACS_GET(IZERO, IZERO, ICONTXT)

  ! Rectangular Grid,  P x 1

  CALL BLACS_GRIDINIT(ICONTXT, ORDER, NP, IONE)
  CALL BLACS_GRIDINFO(ICONTXT, NPROW, NPCOL, MYPROW, MYPCOL)

  !
  !  Get parameters
  !
  CALL GET_PARMS(ICONTXT,CMETHD,PREC,AFMT,IDIM,ISTOPC,ITMAX,ITRACE,ML)
  
  !
  !  Allocate and fill in the coefficient matrix, RHS and initial guess 
  !

  CALL BLACS_BARRIER(ICONTXT,'All')
  T1 = MPI_WTIME()
  CALL CREATE_MATRIX(IDIM,A,B,X,DESC_A,PART_BLOCK,ICONTXT,AFMT)  
  T2 = MPI_WTIME() - T1

  DIM=SIZE(A%ASPK)

  ALLOCATE(H%ASPK(DIM),H%IA1(DIM),H%IA2(DIM),H%PL(SIZE(A%PL)),&
       & H%PL(SIZE(A%PL)),D(SIZE(A%PL)),&
       & DESC_A_OUT%MATRIX_DATA(SIZE(DESC_A%MATRIX_DATA)),&
       & DESC_A_OUT%HALO_INDEX(SIZE(DESC_A%HALO_INDEX)),&
       & DESC_A_OUT%OVRLAP_INDEX(SIZE(DESC_A%OVRLAP_INDEX)),&
       & DESC_A_OUT%OVRLAP_ELEM(SIZE(DESC_A%OVRLAP_ELEM)),&
       & DESC_A_OUT%LOC_TO_GLOB(SIZE(DESC_A%LOC_TO_GLOB)),&
       & DESC_A_OUT%GLOB_TO_LOC(SIZE(DESC_A%GLOB_TO_LOC)), WORK(1024))
  check_descr=15
!  work(5)=9
!!$  WRITE(0,*)'CALLING VERIFY'
!!$  CALL F90_PSVERIFY(D,A,DESC_A,CHECK_DESCR,CONVERT_DESCR,H,&
!!$       & DESC_A_OUT,WORK)
!!$  WRITE(0,*)'VERIFY DONE',CONVERT_DESCR

  deallocate(work)

  CALL DGAMX2D(ICONTXT,'A',' ',IONE, IONE,T2,IONE,T1,T1,-1,-1,-1)
  IF (IAM.EQ.0) Write(6,*) 'Matrix creation Time : ',T2


  !
  !  Prepare the preconditioner.
  !  
  SELECT CASE (PREC)     
  CASE ('SCHW')
     IPREC = 3
  CASE ('ILU')
     IPREC = 2
  CASE ('DIAGSC')
     IPREC = 1
  CASE ('NONE')
     IPREC = 0
  CASE DEFAULT
     WRITE(0,*) 'Unknown preconditioner'
     CALL BLACS_ABORT(ICONTXT,-1)
  END SELECT
  pre%prec=iprec
  CALL BLACS_BARRIER(ICONTXT,'All')
  T1 = MPI_WTIME()
  CALL PRECONDITIONER(A,PRE,DESC_A,IRET)
!!$  CALL PRECONDITIONER(IPREC,A,L,U,D,DESC_A,IRET)
  TPREC = MPI_WTIME()-T1
  
  CALL DGAMX2D(icontxt,'A',' ',IONE, IONE,TPREC,IONE,t1,t1,-1,-1,-1)
  
  IF (IAM.EQ.0) WRITE(6,*) 'Preconditioner Time : ',TPREC

  IF (IRET.NE.0) THEN
     WRITE(0,*) 'Error on preconditioner',IRET
     CALL BLACS_ABORT(ICONTXT,-1)
     STOP
  END IF
 
  !
  ! Iterative method parameters 
  !
  write(*,*) 'Calling Iterative method', size(b),ml
  CALL BLACS_BARRIER(ICONTXT,'All')
  T1 = MPI_WTIME()  
  EPS   = 1.D-9
  IF (CMETHD.EQ.'BICGSTAB') THEN 
    CALL  F90_BICGSTAB(A,PRE,B,X,EPS,DESC_A,& 
         & ITMAX,ITER,ERR,IERR,ITRACE)     
!!$  ELSE IF (CMETHD.EQ.'BICG') THEN 
!!$    CALL  F90_BICG(A,PRE,B,X,EPS,DESC_A,& 
!!$         & ITMAX,ITER,ERR,IERR,ITRACE)     
!!$  ELSE  IF (CMETHD.EQ.'CGS') THEN 
!!$    CALL  F90_CGS(A,PRE,B,X,EPS,DESC_A,& 
!!$         & ITMAX,ITER,ERR,IERR,ITRACE)     
!!$  ELSE IF (CMETHD.EQ.'BICGSTABL') THEN 
!!$    CALL  F90_BICGSTABL(A,PRE,B,X,EPS,DESC_A,& 
!!$         & ITMAX,ITER,ERR,IERR,ITRACE,ML)     
  ELSE
    write(0,*) 'Unknown method ',cmethd
  end IF
  
  CALL BLACS_BARRIER(ICONTXT,'All')
  T2 = MPI_WTIME() - T1
  CALL DGAMX2D(ICONTXT,'A',' ',IONE, IONE,T2,IONE,T1,T1,-1,-1,-1)

  IF (IAM.EQ.0) THEN
     WRITE(6,*) 'Time to Solve Matrix : ',T2
     WRITE(6,*) 'Time per iteration : ',T2/ITER
     WRITE(6,*) 'Number of iterations : ',ITER
     WRITE(6,*) 'Error on exit : ',ERR
     WRITE(6,*) 'INFO  on exit : ',IERR
  END IF

  !  
  !  Cleanup storage and exit
  !
  CALL F90_PSDSFREE(B,DESC_A)
  CALL F90_PSDSFREE(X,DESC_A)
!!$  CALL F90_PSDSFREE(D,DESC_A)

  CALL F90_PSSPFREE(A,DESC_A)
!!$  CALL F90_PSSPFREE(L,DESC_A) 
!!$  CALL F90_PSSPFREE(U,DESC_A)          
  CALL F90_PSDSCFREE(DESC_A,info)
  
  CALL BLACS_GRIDEXIT(ICONTXT)
  CALL BLACS_EXIT(0)
  
  STOP
  
CONTAINS
  !
  ! Get iteration parameters from the command line
  !
  SUBROUTINE  GET_PARMS(ICONTXT,CMETHD,PREC,AFMT,IDIM,ISTOPC,ITMAX,ITRACE,ML)
    integer      :: icontxt
    Character    :: CMETHD*10, PREC*10, AFMT*5
    Integer      :: IDIM, IRET, ISTOPC,ITMAX,ITRACE,ML
    Character*40 :: CHARBUF
    INTEGER      :: IARGC, NPROW, NPCOL, MYPROW, MYPCOL
    EXTERNAL     IARGC
    INTEGER      :: INTBUF(10), IP
    
    CALL BLACS_GRIDINFO(ICONTXT, NPROW, NPCOL, MYPROW, MYPCOL)

    IF (MYPROW==0) THEN
       READ(*,*) IP
       IF (IP.GE.3) THEN
          READ(*,*) CMETHD
          READ(*,*) PREC
          READ(*,*) AFMT
         
        ! Convert strings in array
          DO I = 1, LEN(CMETHD)
             INTBUF(I) = IACHAR(CMETHD(I:I))
          END DO
        ! Broadcast parameters to all processors
          CALL IGEBS2D(ICONTXT,'ALL',' ',10,1,INTBUF,10)
        
          DO I = 1, LEN(PREC)
             INTBUF(I) = IACHAR(PREC(I:I))
          END DO
        ! Broadcast parameters to all processors
          CALL IGEBS2D(ICONTXT,'ALL',' ',10,1,INTBUF,10)
        
          DO I = 1, LEN(AFMT)
             INTBUF(I) = IACHAR(AFMT(I:I))
          END DO
        ! Broadcast parameters to all processors
          CALL IGEBS2D(ICONTXT,'ALL',' ',10,1,INTBUF,10)
        
          READ(*,*) IDIM
          IF (IP.GE.4) THEN
             READ(*,*) ISTOPC
          ELSE
             ISTOPC=1        
          ENDIF
          IF (IP.GE.5) THEN
             READ(*,*) ITMAX
          ELSE
             ITMAX=500
          ENDIF
          IF (IP.GE.6) THEN
             READ(*,*) ITRACE
          ELSE
             ITRACE=-1
          ENDIF
          IF (IP.GE.7) THEN
             READ(*,*) ML
          ELSE
             ML=1
          ENDIF
        ! Broadcast parameters to all processors    
          
          INTBUF(1) = IDIM
          INTBUF(2) = ISTOPC
          INTBUF(3) = ITMAX
          INTBUF(4) = ITRACE
          INTBUF(5) = ML
          CALL IGEBS2D(ICONTXT,'ALL',' ',5,1,INTBUF,5)

          WRITE(6,*)'Solving matrix: ELL1'      
          WRITE(6,*)'on  grid',IDIM,'x',IDIM,'x',IDIM
          WRITE(6,*)' with BLOCK data distribution, NP=',Np,&
               & ' Preconditioner=',PREC,&
               & ' Iterative methd=',CMETHD      
       ELSE
        ! Wrong number of parameter, print an error message and exit
          CALL PR_USAGE(0)      
          CALL BLACS_ABORT(ICONTXT,-1)
          STOP 1
       ENDIF
    ELSE
   ! Receive Parameters
       CALL IGEBR2D(ICONTXT,'ALL',' ',10,1,INTBUF,10,0,0)
       DO I = 1, 10
          CMETHD(I:I) = ACHAR(INTBUF(I))
       END DO
       CALL IGEBR2D(ICONTXT,'ALL',' ',10,1,INTBUF,10,0,0)
       DO I = 1, 10
          PREC(I:I) = ACHAR(INTBUF(I))
       END DO
       CALL IGEBR2D(ICONTXT,'ALL',' ',10,1,INTBUF,10,0,0)
       DO I = 1, 5
          AFMT(I:I) = ACHAR(INTBUF(I))
       END DO
       CALL IGEBR2D(ICONTXT,'ALL',' ',5,1,INTBUF,5,0,0)
       IDIM    = INTBUF(1)
       ISTOPC  = INTBUF(2)
       ITMAX   = INTBUF(3)
       ITRACE  = INTBUF(4)
       ML      = INTBUF(5)
    END IF
    RETURN
    
  END SUBROUTINE GET_PARMS
  !
  !  Print an error message 
  !  
  SUBROUTINE PR_USAGE(IOUT)
    INTEGER :: IOUT
    WRITE(IOUT,*)'Incorrect parameter(s) found'
    WRITE(IOUT,*)' Usage:  pde90 methd prec dim &
         &[istop itmax itrace]'  
    WRITE(IOUT,*)' Where:'
    WRITE(IOUT,*)'     methd:    CGSTAB TFQMR CGS' 
    WRITE(IOUT,*)'     prec :    ILU DIAGSC NONE'
    WRITE(IOUT,*)'     dim       number of points along each axis'
    WRITE(IOUT,*)'               the size of the resulting linear '
    WRITE(IOUT,*)'               system is dim**3'
    WRITE(IOUT,*)'     istop     Stopping criterion  1, 2 or 3 [1]  '
    WRITE(IOUT,*)'     itmax     Maximum number of iterations [500] '
    WRITE(IOUT,*)'     itrace    0  (no tracing, default) or '  
    WRITE(IOUT,*)'               >= 0 do tracing every ITRACE'
    WRITE(IOUT,*)'               iterations ' 
  END SUBROUTINE PR_USAGE

!
!  Subroutine to allocate and fill in the coefficient matrix and
!  the RHS. 
!
  SUBROUTINE CREATE_MATRIX(IDIM,A,B,T,DESC_A,PARTS,ICONTXT,AFMT)
    !
    !   Discretize the partial diferential equation
    ! 
    !   b1 dd(u)  b2 dd(u)    a1 d(u)   a2 d(u)  
    ! -   ------ -  ------ -  -----  -  ------   + a4 u 
    !      dxdx     dydy         dx       dy     
    !
    !  = 0 
    ! 
    ! boundary condition: Dirichlet
    !    0< x,y<1
    !  
    !  u(x,y)(2b1+2b2+a1+a2)+u(x-1,y)(-b1-a1)+u(x,y-1)(-b2-a2)+
    !  -u(x+1,y,z)b1-u(x,y+1,z)b2

    USE TYPESP
    USE TYPEDESC
    USE F90TOOLS
    Implicit None
    INTEGER                  :: IDIM
!!$    external               parts
    integer, parameter       :: nbmax=10
    Real(Kind(1.D0)),Pointer :: B(:),T(:)
    Type (desc_type)  :: DESC_A
    Integer                  :: ICONTXT
    INTERFACE 
      !   .....user passed subroutine.....
      SUBROUTINE PARTS(GLOBAL_INDX,N,NP,PV,NV)
        INTEGER, INTENT(IN)  :: GLOBAL_INDX, N, NP
        INTEGER, INTENT(OUT) :: NV, PV(*) 
      END SUBROUTINE PARTS
    END INTERFACE   ! Local variables
    Type(D_SPMAT)            :: A
    Real(Kind(1.d0))         :: ZT(NBMAX),GLOB_X,GLOB_Y,GLOB_Z
    Integer                  :: M,N,NNZ,GLOB_ROW,J
    Type (D_SPMAT)           :: ROW_MAT
    Integer                  :: X,Y,Z,COUNTER,IA,I,INDX_OWNER
    INTEGER                  :: NPROW,NPCOL,MYPROW,MYPCOL
    Integer                  :: ELEMENT
    INTEGER                  :: INFO, NV, INV
    INTEGER, ALLOCATABLE     :: PRV(:)
    INTEGER, pointer         :: ierrv(:)
    Real(Kind(1.d0)), pointer ::  DWORK(:)
    INTEGER,POINTER        ::  IWORK(:)
    character              :: afmt*5
    ! deltah dimension of each grid cell
    ! deltat discretization time
    Real(Kind(1.D0))         :: DELTAH
    Real(Kind(1.d0)),Parameter   :: RHS=0.d0,ONE=1.d0,ZERO=0.d0
    Real(Kind(1.d0))   :: MPI_WTIME, T1, T2, T3, TINS
    Real(Kind(1.d0))   :: a1, a2, a3, a4, b1, b2, b3 
    external            mpi_wtime,a1, a2, a3, a4, b1, b2, b3
    integer            :: nb, ir1, ir2, ipr
    logical            :: own
    ! common area


    CALL BLACS_GRIDINFO(ICONTXT, NPROW, NPCOL, MYPROW, MYPCOL)

    DELTAH = 1.D0/(IDIM-1)

    ! Initialize array descriptor and sparse matrix storage. Provide an
    ! estimate of the number of non zeroes 
    CALL SETERR(2)
    allocate(ierrv(6))

    ierrv(:) = 0
    M   = IDIM*IDIM
    N   = M
    NNZ = ((N*6)/(NPROW*NPCOL))
    write(*,*) 'Size: n ',n
    Call F90_PSDSCALL(N,N,PARTS,ICONTXT,IERRV,DESC_A)
    write(*,*) 'Allocating A : nnz',nnz
    Call F90_PSSPALL(A,IERRV,DESC_A,NNZ=NNZ)
    ! Define  RHS from boundary conditions; also build initial guess 
    write(*,*) 'Allocating B'
    Call F90_PSDSALL(N,B,IERRV,DESC_A)
    write(*,*) 'Allocating T'
    Call F90_PSDSALL(N,T,IERRV,DESC_A)

    ! We build an auxiliary matrix consisting of one row at a
    ! time; just a small matrix. Might be extended to generate 
    ! a bunch of rows per call. 
    ! 
    ROW_MAT%DESCRA(1:1) = 'G'
    ROW_MAT%FIDA        = 'CSR'
    write(*,*) 'Allocating ROW_MAT',20*nbmax
    ALLOCATE(ROW_MAT%ASPK(20*nbmax),ROW_MAT%IA1(20*nbmax),&
         &ROW_MAT%IA2(20*nbmax),PRV(NPROW),stat=info)
    if (info.ne.0 ) then 
      write(*,*) 'Memory allocation error'
      call blacs_abort(icontxt,-1)
    endif

    TINS = 0.D0
    CALL BLACS_BARRIER(ICONTXT,'ALL')
    T1 = MPI_WTIME()

    ! Loop over rows belonging to current process in a BLOCK
    ! distribution.

    Z=0
    ROW_MAT%IA2(1)=1    
    DO GLOB_ROW = 1, N
      CALL PARTS(GLOB_ROW,N,NPROW,PRV,NV)
      DO INV = 1, NV
        INDX_OWNER = PRV(INV)
        IF (INDX_OWNER == MYPROW) THEN
          ! Local matrix pointer 
          ELEMENT=1
          ! Compute gridpoint Coordinates
          IF (MOD(GLOB_ROW,(IDIM)).EQ.0) THEN
            X = GLOB_ROW/(IDIM)
          ELSE
            X = GLOB_ROW/(IDIM)+1
          ENDIF
          Y = GLOB_ROW-(X-1)*IDIM
          ! GLOB_X, GLOB_Y, GLOB_X coordinates
          GLOB_X=X*DELTAH
          GLOB_Y=Y*DELTAH
          GLOB_Z=Z*DELTAH


          ! Check on boundary points 
          IF (X.EQ.1) THEN
            ROW_MAT%ASPK(ELEMENT)=ONE
            ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM+(Y)
            ELEMENT=ELEMENT+1
          ELSE IF (Y.EQ.1) THEN
            ROW_MAT%ASPK(ELEMENT)=ONE
            ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM+(Y)
            ELEMENT=ELEMENT+1
          ELSE IF (X.EQ.IDIM) THEN
            ROW_MAT%ASPK(ELEMENT)=ONE
            ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM+(Y)
            ELEMENT=ELEMENT+1
          ELSE IF (Y.EQ.IDIM) THEN
            ROW_MAT%ASPK(ELEMENT)=ONE
            ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM+(Y)
            ELEMENT=ELEMENT+1
          ELSE
            ! Internal point: build discretization
            !   
            !  Term depending on   (x-1,y)
            !
            ROW_MAT%ASPK(ELEMENT)=-B1(GLOB_X,GLOB_Y,GLOB_Z)&
                 & -A1(GLOB_X,GLOB_Y,GLOB_Z)
            ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*&
                 & DELTAH)
            ROW_MAT%IA1(ELEMENT)=(X-2)*IDIM+(Y)
            ELEMENT=ELEMENT+1
            !  Term depending on     (x,y-1,z)
            ROW_MAT%ASPK(ELEMENT)=-B2(GLOB_X,GLOB_Y,GLOB_Z)&
                 & -A2(GLOB_X,GLOB_Y,GLOB_Z)
            ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*&
                 & DELTAH)
            ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM+(Y-1)
            ELEMENT=ELEMENT+1
!!$            !  Term depending on     (x,y,z-1)
!!$            ROW_MAT%ASPK(ELEMENT)=-B3(GLOB_X,GLOB_Y,GLOB_Z)&
!!$                 & -A3(GLOB_X,GLOB_Y,GLOB_Z)
!!$            ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*&
!!$                 & DELTAH)
!!$            ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z-1)
!!$            ELEMENT=ELEMENT+1
            !  Term depending on     (x,y,z)
            ROW_MAT%ASPK(ELEMENT)=2*B1(GLOB_X,GLOB_Y,GLOB_Z)&
                 & +2*B2(GLOB_X,GLOB_Y,GLOB_Z)&
                 & +A1(GLOB_X,GLOB_Y,GLOB_Z)&
                 & +A2(GLOB_X,GLOB_Y,GLOB_Z)
            ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*&
                 & DELTAH)
            ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM+(Y)
            ELEMENT=ELEMENT+1                  
!!$            !  Term depending on     (x,y,z+1)
!!$            ROW_MAT%ASPK(ELEMENT)=-B1(GLOB_X,GLOB_Y,GLOB_Z)
!!$            ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*&
!!$                 & DELTAH)
!!$            ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z+1)
!!$            ELEMENT=ELEMENT+1
!!$            !  Term depending on     (x,y+1,z)
            ROW_MAT%ASPK(ELEMENT)=-B2(GLOB_X,GLOB_Y,GLOB_Z)
            ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*&
                 & DELTAH)
            ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM+(Y+1)
            ELEMENT=ELEMENT+1
            !  Term depending on     (x+1,y,z)
            ROW_MAT%ASPK(ELEMENT)=-B3(GLOB_X,GLOB_Y,GLOB_Z)
            ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*&
                 & DELTAH)
            ROW_MAT%IA1(ELEMENT)=(X)*IDIM+(Y)
            ELEMENT=ELEMENT+1
          ENDIF
          ROW_MAT%M=1
          ROW_MAT%K=N
          ROW_MAT%IA2(2)=ELEMENT       
          ! IA== GLOBAL ROW INDEX
          IA=GLOB_ROW
!!$             IA=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z)
!!$        write(0,*) 'Inserting row ',ia,' On proc',myprow
          T3 = MPI_WTIME()
          CALL F90_PSSPINS(A,IA,1,ROW_MAT,IERRV,DESC_A)       
          if (ierrv(1).ne.0) then 
            write(0,*) 'On row ',ia,' IERRV:',ierrv(:)
          endif
          TINS = TINS + (MPI_WTIME()-T3)
          ! Build RHS  
          IF (X==1) THEN     
            GLOB_Y=(Y-IDIM/2)*DELTAH
            GLOB_Z=(Z-IDIM/2)*DELTAH        
            ZT(1) = EXP(-GLOB_Y**2-GLOB_Z**2)
          ELSE IF ((Y==1).OR.(Y==IDIM).OR.(Z==1).OR.(Z==IDIM)) THEN 
            GLOB_X=3*(X-1)*DELTAH
            GLOB_Y=(Y-IDIM/2)*DELTAH
            GLOB_Z=(Z-IDIM/2)*DELTAH        
            ZT(1) = EXP(-GLOB_Y**2-GLOB_Z**2)*EXP(-GLOB_X)
          ELSE
            ZT(1) = 0.D0
          ENDIF
          CALL F90_PSDSINS(1,B,IA,ZT(1:1),IERRV,DESC_A)
          ZT(1)=0.D0
          CALL F90_PSDSINS(1,T,IA,ZT(1:1),IERRV,DESC_A)
        END IF
      END DO
    END DO

    CALL BLACS_BARRIER(ICONTXT,'ALL')    
    T2 = MPI_WTIME()

    WRITE(*,*) '     pspins  time',TINS
    WRITE(*,*) '   Insert time',(T2-T1)

    DEALLOCATE(ROW_MAT%ASPK,ROW_MAT%IA1,ROW_MAT%IA2)

    write(*,*) 'Calling SPASB'
    CALL BLACS_BARRIER(ICONTXT,'ALL')
    T1 = MPI_WTIME()

    CALL F90_PSSPASB(A,IERRV,DESC_A,AFMT=AFMT) 

    CALL BLACS_BARRIER(ICONTXT,'ALL')
    T2 = MPI_WTIME()

    WRITE(0,*) '   Assembly  time',(T2-T1),' ',a%fida(1:4)

    CALL F90_PSDSASB(B,IERRV,DESC_A)
    CALL F90_PSDSASB(T,IERRV,DESC_A)
    IF (MYPROW.EQ.0) THEN
      WRITE(0,*) '   End CREATE_MATRIX'
    ENDIF
    RETURN

  END SUBROUTINE CREATE_MATRIX
END PROGRAM PP2D
!
! Functions parametrizing the differential equation 
!  
FUNCTION A1(X,Y,Z)
  REAL(KIND(1.D0)) :: A1
  REAL(KIND(1.D0)) :: X,Y,Z
  A1=1.D0
END FUNCTION A1
FUNCTION A2(X,Y,Z)
  REAL(KIND(1.D0)) ::  A2
  REAL(KIND(1.D0)) :: X,Y,Z
  A2=2.D1*Y
END FUNCTION A2
FUNCTION A3(X,Y,Z)
  REAL(KIND(1.D0)) ::  A3
  REAL(KIND(1.D0)) :: X,Y,Z      
  A3=1.D0
END FUNCTION A3
FUNCTION A4(X,Y,Z)
  REAL(KIND(1.D0)) ::  A4
  REAL(KIND(1.D0)) :: X,Y,Z      
  A4=1.D0
END FUNCTION A4
FUNCTION B1(X,Y,Z)
  REAL(KIND(1.D0)) ::  B1   
  REAL(KIND(1.D0)) :: X,Y,Z
  B1=1.D0
END FUNCTION B1
FUNCTION B2(X,Y,Z)
  REAL(KIND(1.D0)) ::  B2
  REAL(KIND(1.D0)) :: X,Y,Z
  B2=1.D0
END FUNCTION B2
FUNCTION B3(X,Y,Z)
  REAL(KIND(1.D0)) ::  B3
  REAL(KIND(1.D0)) :: X,Y,Z
  B3=1.D0
END FUNCTION B3