* PACKAGE LSVDSL  !" 線型解析 : SVD 法による特異連立方程式の解
*
*"  [HIS] 93/08/12(takepiro) NUMERICAL RECIPIES 使用
*"        93/09/27(takepiro) 0 固有値判定基準の設定
*
**********************************************************************
      SUBROUTINE LSVDSL
     M            ( AMAT   , 
     M              BVEC   ,
     O              IERR   , 
     I              EPS    ,
     I              MATDIM , OSVDC  )
*
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 
      INCLUDE   (ZLDIM)                      !" 格子点数, 波数
#else
#include        "zcdim.F"                    !" 
#include        "zldim.F"                    !" 係数行列の大きさ
#endif
*
*   [MODIFY] 
      REAL      AMAT  ( MATDMX, MATDMX )     !" 係数行列
      REAL      BVEC  ( MATDMX )             !" 右辺ベクトル / 解
*
*   [OUTPUT] 
      INTEGER   IERR                         !" サブルーチンエラーコード
*
*   [INPUT] 
      REAL      EPS                          !" 特異性判定基準
      INTEGER   MATDIM                       !" 計算する次元
      LOGICAL   OSVDC                        !" SVD 分解計算スイッチ
*
*   [INTERNAL WORK] 
      REAL      XVEC  ( MATDMX )             !" 固有値実数部
      REAL      WDIAG ( MATDMX )             !" 固有値虚数部
      REAL      WORKV ( MATDMX, MATDMX )     !" 固有ベクトル実数部
*
      INTEGER   IUNIT , I
      LOGICAL   OFIRST 
      DATA      OFIRST  / .TRUE. /
*
      SAVE      OFIRST
*
      IF ( OFIRST )THEN
         OFIRST = .FALSE.
         CALL GLPGET( 'MSGUNIT', IUNIT )
         WRITE ( IUNIT, * ) ' SINGULAR VALUE DECOMPOSITION : 93/09/27'
      ENDIF
*
*"  < 1. SVD 分解 >
*
      IF ( OSVDC )THEN
         CALL SVDCMP
     M        ( AMAT   , 
     O          IERR   ,
     I          MATDMX , MATDIM , 
     I          MATDMX , MATDIM , 
     W          WDIAG  , WORKV    )
      ENDIF
*
*"  < 2. 特異性の除去 >
*
      DO 100 I = 1, MATDIM
         IF ( ABS(WDIAG(I)) .LT. EPS )THEN
            WDIAG(I) = 0.
         ENDIF
 100  CONTINUE 
*
*"  < 3. 解を求める >
*
      CALL SVBKSB
     I       ( AMAT   , 
     W         WDIAG  , WORKV  ,
     I         MATDMX , MATDIM , 
     I         MATDMX , MATDIM , 
     I         BVEC   , 
     O         XVEC               )
*
      CALL COPY( BVEC , XVEC , MATDIM )
*
      RETURN
      END