RFLU_InvertMatrixSVD.F90

Go to the documentation of this file.

! *********************************************************************
! * Rocstar Simulation Suite                                          *
! * Copyright@2015, Illinois Rocstar LLC. All rights reserved.        *
! *                                                                   *
! * Illinois Rocstar LLC                                              *
! * Champaign, IL                                                     *
! * www.illinoisrocstar.com                                           *
! * sales@illinoisrocstar.com                                         *
! *                                                                   *
! * License: See LICENSE file in top level of distribution package or *
! * http://opensource.org/licenses/NCSA                               *
! *********************************************************************
! *********************************************************************
! * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,   *
! * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   *
! * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          *
! * NONINFRINGEMENT.  IN NO EVENT SHALL THE CONTRIBUTORS OR           *
! * COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER       *
! * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,   *
! * Arising FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE    *
! * USE OR OTHER DEALINGS WITH THE SOFTWARE.                          *
! *********************************************************************
! ******************************************************************************
!
! Purpose: Invert general m-by-n matrix A using SVD
!
! Description: None.
!
! Input:
!   global      Pointer to global data
!   nRows       Number of rows of A
!   nCols       Number of columns of A
!   a           Matrix A, to be inverted
!
! Output:
!   aInv        Inverse of A   
!   sCount      Number of singular values below threshold 
!
! Notes:
!   1. Uses LAPACK routines to carry out SVD.
!
! ******************************************************************************
!
! $Id: RFLU_InvertMatrixSVD.F90,v 1.9 2008/12/06 08:44:12 mtcampbe Exp $
!
! Copyright: (c) 2002-2005 by the University of Illinois
!
! ******************************************************************************

SUBROUTINE rflu_invertmatrixsvd(global,nRows,nCols,a,aInv,sCount)

  USE moddatatypes
  USE moderror
  USE modglobal, ONLY: t_global
 
  USE modsortsearch
 
  IMPLICIT NONE

! ******************************************************************************
! Definitions and declarations
! ******************************************************************************

! ==============================================================================
! Arguments
! ==============================================================================

  INTEGER, INTENT(IN) :: nrows,ncols
  INTEGER, INTENT(OUT), OPTIONAL :: scount
  REAL(RFREAL) :: a(nrows,ncols),ainv(ncols,nrows)
  TYPE(t_global), POINTER :: global

! ==============================================================================
! Locals
! ==============================================================================

  CHARACTER(CHRLEN) :: rcsidentstring
  INTEGER :: errorflag,i,info,j,workarrayrealsize
  INTEGER :: workarrayint(8*ncols)  
  REAL(RFREAL) :: s(ncols)
  REAL(RFREAL), DIMENSION(:), ALLOCATABLE :: workarrayreal
  REAL(RFREAL) :: sinv(ncols,nrows),u(nrows,nrows),vt(ncols,ncols), & 
                  v(ncols,ncols)

! ******************************************************************************
! Start
! ******************************************************************************

  rcsidentstring = '$RCSfile: RFLU_InvertMatrixSVD.F90,v $ $Revision: 1.9 $'

  CALL registerfunction(global,'RFLU_InvertMatrixSVD',&
  'RFLU_InvertMatrixSVD.F90')

! ******************************************************************************
! Set work array size and allocate memory
! ******************************************************************************

  workarrayrealsize = 2*(4*ncols*ncols + nrows + 9*ncols)

  ALLOCATE(workarrayreal(workarrayrealsize),stat=errorflag)
  global%error = errorflag    
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_allocate,__line__,'workArrayReal')
  END IF ! global%error  

! ******************************************************************************
! Perform SVD
! ******************************************************************************
            
  CALL dgesdd('A',nrows,ncols,a,nrows,s,u,nrows,vt,ncols,workarrayreal, & 
              workarrayrealsize,workarrayint,info)
  global%error = info
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_lapack_output,__line__)
  END IF ! global%error

! ******************************************************************************
! Construct actual singular-value matrix
! ******************************************************************************

  scount = 0
  
  DO i = 1,ncols ! Explicit loop to avoid Frost problem
    DO j = 1,nrows
      sinv(i,j) = 0.0_rfreal
    END DO ! j
  END DO ! i

! Based on absolute size
!  DO i = 1,nCols  
!    IF ( s(i) > 100.0_RFREAL*EPSILON(1.0_RFREAL) ) THEN   
!      sInv(i,i) = 1.0_RFREAL/s(i)
!    ELSE 
!      sInv(i,i) = 0.0_RFREAL
!      sCount    = sCount + 1
!    END IF ! s
!  END DO ! i

! Based on ratio of successive singular values
  IF ( s(1) > 100.0_rfreal*epsilon(1.0_rfreal) ) THEN 
    sinv(1,1) = 1.0_rfreal/s(1)
  ELSE 
    sinv(1,1) = 0.0_rfreal
    scount = 1
  END IF ! s

  outerloop: DO i = 2,ncols  
    IF ( s(i) > 0.01_rfreal*s(i-1) ) THEN   
      sinv(i,i) = 1.0_rfreal/s(i)
    ELSE 
      sinv(i,i) = 0.0_rfreal
      scount    = scount + 1
      
      DO j = i+1,ncols
        sinv(j,j) = 0.0_rfreal
        scount = scount + 1
      END DO ! j
      
      EXIT outerloop
    END IF ! s
  END DO outerloop
 
! ******************************************************************************
! Compute inverse
! ******************************************************************************

  ainv = matmul(transpose(vt),matmul(sinv,transpose(u))) ! Lapack

! ******************************************************************************
! End
! ******************************************************************************

  DEALLOCATE(workarrayreal,stat=errorflag)
  global%error = errorflag    
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_deallocate,__line__,'workArrayReal')
  END IF ! global%error  

  CALL deregisterfunction(global)

END SUBROUTINE rflu_invertmatrixsvd

! ******************************************************************************
!
! RCS Revision history:
!
! $Log: RFLU_InvertMatrixSVD.F90,v $
! Revision 1.9  2008/12/06 08:44:12  mtcampbe
! Updated license.
!
! Revision 1.8  2008/11/19 22:17:25  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.7  2006/04/07 15:19:16  haselbac
! Removed tabs
!
! Revision 1.6  2005/10/05 13:51:42  haselbac
! Cosmetics
!
! Revision 1.5  2003/12/04 03:23:56  haselbac
! Fixed bug, changed init, added new method of detecting singularity
!
! Revision 1.4  2003/07/22 01:56:50  haselbac
! Change to init of sInv and cosmetics
!
! Revision 1.3  2002/10/08 15:48:56  haselbac
! {IO}STAT=global%error replaced by {IO}STAT=errorFlag - SGI problem
!
! Revision 1.2  2002/09/09 14:15:01  haselbac
! global now under regions, bug fix: workArrayReal not deallocated
!
! Revision 1.1  2002/07/25 14:34:59  haselbac
! Initial revision
!
! ******************************************************************************