NSSOLVER.f

!     This software is Copyright (c) 2012-2015 The Regents of the 
!     University of California. All Rights Reserved.
!
!     Permission to copy and modify this software and its documentation
!     for educational, research and non-profit purposes, without fee, 
!     and without a written agreement is hereby granted, provided that
!     the above copyright notice, this paragraph and the following three
!     paragraphs appear in all copies.
!
!     Permission to make commercial use of this software may be obtained
!     by contacting:
!
!     Technology Transfer Office
!     9500 Gilman Drive, Mail Code 0910
!     University of California
!     La Jolla, CA 92093-0910
!     (858) 534-5815
!     invent@ucsd.edu
!
!     This software program and documentation are copyrighted by The
!     Regents of the University of California. The software program and
!     documentation are supplied "as is", without any accompanying
!     services from The Regents. The Regents does not warrant that the
!     operation of the program will be uninterrupted or error-free. The
!     end-user understands that the program was developed for research
!     purposes and is advised not to rely exclusively on the program for
!     any reason.
!
!     IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY 
!     PARTY FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL 
!     DAMAGES, INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS 
!     SOFTWARE AND ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF 
!     CALIFORNIA HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 
!     THE UNIVERSITY OF CALIFORNIA SPECIFICALLY DISCLAIMS ANY 
!     WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 
!     OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE 
!     SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND THE 
!     UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE 
!     MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
      
      SUBROUTINE NSSOLVER(nFaces, gnNo, dof, nNo, nnz, mynNo, commu, 
     2   cS, face, ls, rowPtr, colPtr, Val, Ri)
      
      INCLUDE "svLS_STD.h"
      

      INTEGER, INTENT(IN) :: nFaces, gnNo, dof, nNo, nnz, mynNo
      TYPE(svLS_commuType), INTENT(IN) :: commu
      TYPE(svLS_cSType), INTENT(IN) :: cS(commu%nTasks)
      TYPE(svLS_faceType), INTENT(INOUT) :: face(nFaces)
      TYPE(svLS_lsType), INTENT(INOUT) :: ls
      INTEGER, INTENT(IN) :: rowPtr(2,nNo), colPtr(nnz)
      REAL*8, INTENT(IN) :: Val(dof*dof,nnz)
      REAL*8, INTENT(INOUT) :: Ri(dof,nNo)

      INTEGER i, j, k, iB, iBB, nB, nsd
      REAL*8 CPUT, NORMS, NORMV, DOTS, DOTV, eps
      REAL*8, ALLOCATABLE :: U(:,:,:), P(:,:),
     2   MU(:,:,:), MP(:,:), A(:,:), B(:), xB(:), mK(:,:), mG(:,:), 
     3   mD(:,:), mL(:), Gt(:,:), Rm(:,:), Rc(:), Rmi(:,:), Rci(:)
 
      nsd = dof - 1
      iB = ls%RI%mItr
      nB = 2*iB
      ALLOCATE(Rm(nsd,nNo), Rc(nNo), Rmi(nsd,nNo), Rci(nNo), 
     2   U(nsd,nNo,iB), P(nNo,iB), MU(nsd,nNo,nB), MP(nNo,nB), 
     3   A(nB,nB), B(nB), xB(nB))

      Rmi = Ri(1:nsd,:)
      Rci = Ri(dof,:)

      xB          = 0D0
      B           = 0D0
      Rm          = Rmi
      Rc          = Rci
      eps         = SQRT(NORMV(nsd, mynNo, commu, Rm)**2D0
     2            +      NORMS(     mynNo, commu, Rc)**2D0)
      ls%RI%iNorm = eps
      ls%RI%fNorm = eps
      ls%CG%callD = 0D0
      ls%GM%callD = 0D0
      ls%CG%itr   = 0
      ls%GM%itr   = 0
      ls%RI%callD = CPUT()
      ls%RI%suc   = .FALSE.
      eps         = MAX(ls%RI%absTol,ls%RI%relTol*eps)

      CALL DEPART
      CALL BCPRE

      DO i=1, ls%RI%mItr
         iB  = 2*i - 1
         iBB = 2*i
         ls%RI%dB = ls%RI%fNorm
         CALL GMRES(nFaces, nsd, nNo, nnz, mynNo, commu, cS, face, 
     2      ls%GM, rowPtr, colPtr, mK, Rm, U(:,:,i))  
         
         CALL SPARMULVS(nsd, nNo, nnz, commu, cS, 
     2      rowPtr, colPtr, mD, U(:,:,i), P(:,i))   
         
         P(:,i) = Rc - P(:,i)       
         CALL CGRAD(nFaces, nsd, nNo, nnz, mynNo, commu, cS, face, 
     2      ls%CG, rowPtr, colPtr, Gt, mG, mL, P(:,i)) 
         
         CALL SPARMULSV(nsd, nNo, nnz, commu, cS, 
     2      rowPtr, colPtr, mG, P(:,i), MU(:,:,iB)) 
         
         MU(:,:,iBB) = Rm - MU(:,:,iB)   
         CALL GMRES(nFaces, nsd, nNo, nnz, mynNo, commu, cS, face, 
     2      ls%GM, rowPtr, colPtr, mK, MU(:,:,iBB), U(:,:,i)) 
         
         CALL SPARMULVV(nsd, nNo, nnz, commu, cS, 
     2      rowPtr, colPtr, mK, U(:,:,i), MU(:,:,iBB))
         
         CALL ADDBCMUL(BCOP_TYPE_ADD, nFaces, nsd, nNo, mynNo, commu, 
     2      face, U(:,:,i), MU(:,:,iBB))
         
         CALL SPARMULSS(nNo, nnz, commu, cS, 
     2      rowPtr, colPtr, mL, P(:,i), MP(:,iB))
         
         CALL SPARMULVS(nsd, nNo, nnz, commu, cS, 
     2      rowPtr, colPtr, mD, U(:,:,i), MP(:,iBB))

         DO k=iB, iBB
            DO j=1, k - 1
               A(j,k) = DOTV(nsd, mynNo, commu, MU(:,:,j), MU(:,:,k))
     2                + DOTS(     mynNo, commu, MP(:,j),   MP(:,k))
               A(k,j) = A(j,k)
            END DO
            A(k,k) = NORMV(nsd, mynNo, commu, MU(:,:,k))**2D0 
     2             + NORMS(     mynNo, commu, MP(:,k))**2D0
            B(k)   = DOTV (nsd, mynNo, commu, MU(:,:,k), Rmi) 
     2             + DOTS (     mynNo, commu, MP(:,k), Rci)
         END DO

         xB = B
         CALL GE(iBB, A(1:iBB,1:iBB), xB(1:iBB))

         ls%RI%fNorm = SQRT(ls%RI%iNorm**2D0 - SUM(xB(1:iBB)*B(1:iBB)))
         IF(ls%RI%fNorm .LT. eps) THEN
            ls%RI%suc = .TRUE.
            EXIT
         END IF
 
         Rm = Rmi - xB(1)*MU(:,:,1)
         Rc = Rci - xB(1)*MP(:,1)
         DO j=2, iBB
            Rm = Rm - xB(j)*MU(:,:,j)
            Rc = Rc - xB(j)*MP(:,j)
         END DO
      END DO
      IF (i .GT. ls%RI%mItr) THEN
         ls%RI%itr = ls%RI%mItr
      ELSE
         ls%RI%itr = i
 
         Rc = Rci - xB(1)*MP(:,1)
         DO j=2, iBB
            Rc = Rc - xB(j)*MP(:,j)
         END DO
      END IF
      ls%Resc = NINT(1D2*(NORMS(mynNo, commu, Rc)/
     2   ls%RI%fNorm)**2D0)
      ls%Resm = 100 - ls%Resc
 
      Rmi = xB(2)*U(:,:,1)
      Rci = xB(1)*P(:,1)
      DO i=2, ls%RI%itr
         iB  = 2*i - 1
         iBB = 2*i

         Rmi = Rmi + xB(iBB)*U(:,:,i)
         Rci = Rci + xB(iB)*P(:,i)
      END DO

      ls%RI%callD = CPUT() - ls%RI%callD
      ls%RI%dB    = 1D1*LOG(ls%RI%fNorm/ls%RI%dB)
      
      Ri(1:nsd,:) = Rmi
      Ri(dof,:) = Rci
      
      DEALLOCATE (Rm, Rc, Rmi, Rci, U, P, MU, MP, A, B, mK, mD, mG, mL, 
     2   Gt)
 
      IF (commu%masF) CALL LOGFILE

      RETURN
      CONTAINS

!====================================================================

      SUBROUTINE DEPART

      IMPLICIT NONE

      INTEGER i, j, k, l
      REAL*8, ALLOCATABLE :: tmp(:)
      
      ALLOCATE(mK(nsd*nsd,nnz), mG(nsd,nnz), mD(nsd,nnz), mL(nnz), 
     2   Gt(nsd,nnz), tmp((nsd+1)*(nsd+1)))

      IF (nsd .EQ. 2) THEN
         DO i=1, nnz
            tmp = Val(:,i)

            mK(1,i) = tmp(1)
            mK(2,i) = tmp(2)
            mK(3,i) = tmp(4)
            mK(4,i) = tmp(5)

            mG(1,i) = tmp(3)
            mG(2,i) = tmp(6)

            mD(1,i) = tmp(7)
            mD(2,i) = tmp(8)
            
            mL(i)   = tmp(9)
         END DO
      ELSE IF(nsd .EQ. 3) THEN
         DO i=1, nnz
            tmp = Val(:,i)
   
            mK(1,i) = tmp(1)
            mK(2,i) = tmp(2)
            mK(3,i) = tmp(3)
            mK(4,i) = tmp(5)
            mK(5,i) = tmp(6)
            mK(6,i) = tmp(7)
            mK(7,i) = tmp(9)
            mK(8,i) = tmp(10)
            mK(9,i) = tmp(11)
            
            mG(1,i) = tmp(4)
            mG(2,i) = tmp(8)
            mG(3,i) = tmp(12)

            mD(1,i) = tmp(13)
            mD(2,i) = tmp(14)
            mD(3,i) = tmp(15)

            mL(i)   = tmp(16)
         END DO
      ELSE
         PRINT *, "Not defined nsd for DEPART", nsd
      END IF
 
      DO i=1, nNo
         Do j=rowPtr(1,i), rowPtr(2,i)
            k = colPtr(j)
            DO l=rowPtr(1,k), rowPtr(2,k)
               IF (colPtr(l) .EQ. i) THEN
                  Gt(:,l) = -mG(:,j)
                  EXIT
               END IF
            END DO
         END DO
      END DO

      RETURN
      END SUBROUTINE DEPART

!====================================================================
      
      SUBROUTINE BCPRE

      IMPLICIT NONE

      INTEGER faIn, i, a, Ac
      REAL*8 NORMV
      REAL*8, ALLOCATABLE :: v(:,:)

      DO faIn=1, nFaces
         IF (face(faIn)%coupledFlag) THEN
            IF (face(faIn)%sharedFlag) THEN
               IF (.NOT.ALLOCATED(v)) ALLOCATE(v(nsd,nNo))
               v = 0D0
               DO a=1, face(faIn)%nNo
                  Ac = face(faIn)%glob(a)
                  DO i=1, nsd
                     v(i,Ac) = face(faIn)%valM(i,a)
                  END DO
               END DO
               face(faIn)%nS = NORMV(nsd, mynNo, commu, v)**2D0
            ELSE
               face(faIn)%nS = 0D0
               DO a=1, face(faIn)%nNo
                  Ac = face(faIn)%glob(a)
                  DO i=1, nsd
                     face(faIn)%nS = face(faIn)%nS + 
     2                  face(faIn)%valM(i,a)**2D0
                  END DO
               END DO
            END IF
         END IF
      END DO

      RETURN
      END SUBROUTINE BCPRE

!====================================================================

      SUBROUTINE LOGFILE

      IMPLICIT NONE
 
      LOGICAL flag
      INTEGER fid, i, j
      CHARACTER*16, PARAMETER :: fName = 'svLS_NS.log'

      INQUIRE(FILE=fName, EXIST=flag)

      fid = 11232
      OPEN(fid, FILE=fName, POSITION='APPEND')
      
      IF (.NOT.flag) THEN
         i = 0
         DO j=1, nFaces
            IF (face(j)%coupledFlag) i = i + 1
         END DO
         WRITE(fid,*) gnNo, commu%nTasks, i
      END IF

      i = 0
      IF (ls%RI%suc) i = i + 100
      IF (ls%GM%suc) i = i + 10
      IF (ls%CG%suc) i = i + 1

      WRITE(fid,"(I4.3,I3,I4,I5,3I4,3ES9.2E2,3I4)") 
     2   i, ls%RI%itr, ls%GM%itr, ls%CG%itr, 
     3   NINT((ls%RI%CallD-ls%GM%CallD-ls%CG%CallD)/ls%RI%CallD*1D2),
     4   NINT(ls%GM%callD/ls%RI%CallD*1D2),
     5   NINT(ls%CG%callD/ls%RI%CallD*1D2),
     6   ls%RI%iNorm, ls%RI%fNorm/ls%RI%iNorm, ls%RI%CallD,
     7   ls%Resm, ls%Resc, NINT(ls%RI%dB) 

      CLOSE(fid)

      RETURN
      END SUBROUTINE LOGFILE

      END SUBROUTINE NSSOLVER