Commit e20a3def authored by Nichols, Stephen's avatar Nichols, Stephen
Browse files

Streamlined the routines in matrix_ops.f90 : 

    (a) moved the indexing of the arrays to the subroutine calls in custom_linsolve.f90
    (b) this removed all the "i,j,k" indexing from the routines in matrix_ops.f90
parent 84d8c7a0

custom_linsolve.f90

+2 −2
Original line number Diff line number Diff line
@@ -318,7 +318,7 @@ program custom_linsolve 
  do m = 1,nz

     do l = 1,ny

        do k = 1,nx

           call matrix_invert(n,k,l,m,matrix_info)

           call matrix_invert(n,pivot(:,k,l,m),A(:,:,k,l,m),Ao(:,:,k,l,m),z(:,k,l,m),matrix_info)

        enddo

     enddo

  enddo
@@ -575,7 +575,7 @@ program custom_linsolve 
              ! perform lu decomp and back substitution ... returns z

              !! yes, we've already performed the lu decomp for the matrix

              !! inversion, but I'm doing it again to check for correctness

              call solve(n,k,l,m,matrix_info)

              call solve(n,pivot(:,k,l,m),A(:,:,k,l,m),z(:,k,l,m),matrix_info)



           enddo

        enddo

matrix_ops.f90

+52 −38
Original line number Diff line number Diff line
@@ -263,14 +263,17 @@ 
      end subroutine check_matrix



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

      subroutine fwd_bk_sub(neqn,nn,i,j,k)

      subroutine fwd_bk_sub(neqn,nn,ppivot,aa,zz)

      implicit none



#ifdef USE_OPENMP_OFFLOAD

      !$omp declare target

#endif



      integer, intent(in) :: neqn,nn,i,j,k

      integer, intent(in) :: neqn,nn

      integer, intent(inout) :: ppivot(neqn)

      real *8, intent(inout) :: aa(neqn,neqn)

      real *8, intent(inout) :: zz(neqn)

      real*8 :: zpl,zpm

      integer :: l,m

      integer :: pl,pm
@@ -278,13 +281,13 @@ 
      ! perform the same eliminations on Z that were

      ! completed on A (this is the quasi-forward substitution pass)

      do m = nn,neqn-1

          pm = pivot(m,i,j,k)

          if (z(pm,i,j,k) .ne. 0.d0) then    !! no need to mult by 0.0

              zpm = z(pm,i,j,k)

          pm = ppivot(m)

          if (zz(pm) .ne. 0.d0) then    !! no need to mult by 0.0

              zpm = zz(pm)

              do l = m+1,neqn

                  pl = pivot(l,i,j,k)

                  if (a(pl,m,i,j,k) .ne. 0.d0) then !! no need to add 0.0

                      z(pl,i,j,k) = z(pl,i,j,k) - a(pl,m,i,j,k)*zpm

                  pl = ppivot(l)

                  if (aa(pl,m) .ne. 0.d0) then !! no need to add 0.0

                      zz(pl) = zz(pl) - aa(pl,m)*zpm

                  endif

              end do

          endif
@@ -293,14 +296,14 @@ 
      ! solves the upper triangular system using column

      ! oriented back substitution

      do l = neqn,1,-1

          pl = pivot(l,i,j,k)

          if (z(pl,i,j,k) .ne. 0.d0) then     !! no need to mult by 0.0

              z(pl,i,j,k) = z(pl,i,j,k)/a(pl,l,i,j,k)

              zpl = z(pl,i,j,k)

          pl = ppivot(l)

          if (zz(pl) .ne. 0.d0) then     !! no need to mult by 0.0

              zz(pl) = zz(pl)/aa(pl,l)

              zpl = zz(pl)

              do m = l-1,1,-1

                  pm = pivot(m,i,j,k)

                  if (a(pm,l,i,j,k) .ne. 0.d0) then !! no need to add 0.0

                      z(pm,i,j,k) = z(pm,i,j,k) - a(pm,l,i,j,k)*zpl

                  pm = ppivot(m)

                  if (aa(pm,l) .ne. 0.d0) then !! no need to add 0.0

                      zz(pm) = zz(pm) - aa(pm,l)*zpl

                  endif

              end do

          endif
@@ -310,14 +313,18 @@ 
      end subroutine fwd_bk_sub



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

      subroutine find_inverse(neqn,i,j,k)

      subroutine find_inverse(neqn,ppivot,aa,aao,zz)

      implicit none



#ifdef USE_OPENMP_OFFLOAD

      !$omp declare target

#endif



      integer, intent(in) :: neqn,i,j,k

      integer, intent(in) :: neqn

      integer, intent(inout) :: ppivot(neqn)

      real *8, intent(inout) :: aa(neqn,neqn)

      real *8, intent(inout) :: aao(neqn,neqn)

      real *8, intent(inout) :: zz(neqn)

      integer :: l,m

      integer :: pm
@@ -332,30 +339,30 @@ 
      do l = 1,neqn

          ! initialize z with the appropriate column of the Identity matrix

          do m = 1, neqn

              z(m,i,j,k) = 0.d0

              zz(m) = 0.d0

          enddo

          z(pivot(l,i,j,k),i,j,k) = 1.d0

          zz(ppivot(l)) = 1.d0



          ! solve for the appropriate column of the inverse of a

          ! NOTE: uses the pivoted a matrix

          ! NOTE2: can send in "l" iff we're initializing z with 

          !        a column of the Identity matrix ... otherwise

          !        we must send in ONE : "1" instead of "l"

          call fwd_bk_sub(neqn,l,i,j,k)

          call fwd_bk_sub(neqn,l,ppivot,aa,zz)

          

          ! store the column in a

          ! NOTE: pivoting is reversed so that the matrix returned in a

          !       is truly the inverse of the original a

          ! NOTE2: also need to pivot the column with this algorithm

          do m = 1,neqn

              pm = pivot(m,i,j,k)

              ao(m,pivot(l,i,j,k),i,j,k) = z(pm,i,j,k)

              pm = ppivot(m)

              aao(m,ppivot(l)) = zz(pm)

          enddo

      enddo



      do m = 1,neqn

         do l = 1,neqn

            a(l,m,i,j,k) = ao(l,m,i,j,k)

            aa(l,m) = aao(l,m)

         enddo

      enddo
@@ -364,26 +371,30 @@ 


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



      subroutine matrix_invert(neqn,i,j,k,matrix_info)

      subroutine matrix_invert(neqn,ppivot,aa,aao,zz,matrix_info)

      implicit none



#ifdef USE_OPENMP_OFFLOAD

      !$omp declare target

#endif



      integer, intent(in) ::  i,j,k,neqn

      integer, intent(in) :: neqn

      integer, intent(inout) :: ppivot(neqn)

      real *8, intent(inout) :: aa(neqn,neqn)

      real *8, intent(inout) :: aao(neqn,neqn)

      real *8, intent(inout) :: zz(neqn)

      integer, intent(inout) :: matrix_info



      ! routines for inverting the C matrix using custom library



      call ludecomp(neqn,pivot(1,i,j,k),a(1,1,i,j,k),matrix_info)

      call ludecomp(neqn,ppivot,aa,matrix_info)



      if (matrix_info .gt. 0) then

         !! bail if there is problem with the matrix

         return

      else

         !! find the inverse if the matrix passes the checks

         call find_inverse(neqn,i,j,k)

         call find_inverse(neqn,ppivot,aa,aao,zz)

      endif



      return
@@ -391,14 +402,17 @@ 


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



  subroutine solve(neqn,i,j,k,matrix_info)

  subroutine solve(neqn,ppivot,aa,zz,matrix_info)

  implicit none



#ifdef USE_OPENMP_OFFLOAD

  !$omp declare target

#endif



  integer, intent(in) ::  neqn,i,j,k

  integer, intent(in) ::  neqn

  integer, intent(inout) :: ppivot(neqn)

  real *8, intent(inout) :: aa(neqn,neqn)

  real *8, intent(inout) :: zz(neqn)

  integer, intent(inout) :: matrix_info



  real*8 :: z_tmp
@@ -407,32 +421,32 @@ 
  ! routines for inverting the C matrix using custom library



  ! lu decomp of matrix c

  call ludecomp(neqn,pivot(1,i,j,k),a(1,1,i,j,k),matrix_info)

  call ludecomp(neqn,ppivot,aa,matrix_info)



  if (matrix_info .gt. 0) then

     !! bail if there is problem with the matrix

     return

  else

     call fwd_bk_sub(neqn,1,i,j,k)

     call fwd_bk_sub(neqn,1,ppivot,aa,zz)



     ! reverse the pivoting without an additional array

     do l = 1,neqn-1

        if (l .ne. pivot(l,i,j,k)) then

        if (l .ne. ppivot(l)) then



           ! find the row to switch

           l_tmp = 0

           do m = l+1, neqn

              if (pivot(m,i,j,k) .eq. l) l_tmp = m

              if (ppivot(m) .eq. l) l_tmp = m

           end do



           ! switch rows in solution vector

           z_tmp = z(pivot(l,i,j,k),i,j,k)

           z(pivot(l,i,j,k),i,j,k) = z(pivot(l_tmp,i,j,k),i,j,k)

           z(pivot(l_tmp,i,j,k),i,j,k) = z_tmp

           z_tmp = zz(ppivot(l))

           zz(ppivot(l)) = zz(ppivot(l_tmp))

           zz(ppivot(l_tmp)) = z_tmp



           ! update pivot vector

           pivot(l_tmp,i,j,k) = pivot(l,i,j,k)

           pivot(l,i,j,k) = l

           ppivot(l_tmp) = ppivot(l)

           ppivot(l) = l

        end if

     end do