!**********************************************************************
! GSM Solver Ver 1.0
! Generalized Scattering Matrix (GSM) Analysis Code
!
!  Copyright (c) 2007 Takuichi Hirano (Tokyo Institute of Technology).
!  All rights reserved.
!                      E-mail: hira@antenna.ee.titech.ac.jp
!                      http://www-antenna.ee.titech.ac.jp/~hira/
!
! [Keywords]
!   * Generalized Scattering Matrix (GSM)
!   * Scattering Matrix (S-Matrix)
!   * Microwave Circuit
!   * Antenna
!   * High Frequency
!   * Electromagnetic Wave
!   * Distributed Constant Circuit
!   * Connection
!   * Mode
!   * Transmission Line
!   * Waveguide
!   * Microstrip Line
!   * Coaxial Cable
!
! [Language and Compiler]
!   Language: Fortran90
!   Compile Environment: Compaq Visual Fortran Version 6.5
!
! [Usage]
!   gsm [input S-matrix file] [input topology file] [output file]
!
! [Acknowledgment]
!   This program uses LAPACK SVD routine for pseudo inverse.
!**********************************************************************

!**********************************************************************
! Modules
!**********************************************************************

!======================================================================
! Constants
!======================================================================
module mod_consts
   real(8),parameter :: eps=1.0d-100                   ! numerically-small number
   real(8),parameter :: eps2=1.0d-16                   ! numerically-small number

   real(8),parameter :: pi=3.141592653589793d0         ! circle ratio
   complex(8),parameter :: zj=dcmplx(0.0d0,1.0d0)      ! pure imaginary
end module

!======================================================================
! GSM (Generalized Scattering Matrix)
!======================================================================
module mod_gsm
   ! **** Type (block, loc_port) ****
   type block_loc_port_list
      integer :: block     ! Block No.
      integer :: loc_port  ! Port No.
   end type
   ! **** Type (gport, in_out) ****
   type gport_in_out
      integer :: gport     ! Block No.
      integer :: in_out    ! Port No.
   end type

   integer :: n_block                  ! No. of Blocks in GSM
   integer,allocatable :: n_port(:)    ! No. of Ports in Each Block
   integer :: m_port                   ! Max. No. of Ports in Each Block
   integer :: n_output_port            ! No. of Ports for Output
   integer :: n_gsm_total_port         ! No. of Ports in GSM System
   integer :: n_gsm_matrix_unknown     ! No. of Unknowns in GSM Matrix

   ! Filename
   character(len=256) :: file_s_matrix,file_topology,file_output

   ! Tables
   type (block_loc_port_list),allocatable :: tab_gport_2_blk_lport(:)
   integer,allocatable :: tab_blk_lport_2_gport(:,:)
   integer,allocatable :: tab_port_connect_info(:)      ! gport -> gport (0: load, -1: excitation)
   complex(8),allocatable :: tab_ex_val_lis(:)          ! Excitation
   type (gport_in_out),allocatable :: tab_out_lis(:)    ! List of Output (c.f. n_output_port)
   integer,allocatable :: tab_unknown_no(:,:)           ! (gport, 1: input; 2: output) -> Unknown No. (0: load, known[=0], -1: excitation, known[=tab_ex_val_lis])

   ! Matrices
   complex(8),allocatable :: s_matrix(:,:,:)    ! S Matrix in Each Block

   complex(8),allocatable :: gsm_matrix(:,:)    ! GSM System Matrix (L.H.S)
   complex(8),allocatable :: gsm_rhs_vec(:)     ! GSM System Vector (R.H.S)
end module

!**********************************************************************
! Main Program
!**********************************************************************
program main
   use dflib    ! built-in module for command-line argument: getarg
   use mod_consts
   use mod_gsm
   implicit none

   integer(2) :: status
   integer :: ierr

   if(nargs()/=4) then
      ! No. of arguments is not correct.
      call show_help
      stop
   end if
   call getarg(1,file_s_matrix,status)
   call getarg(2,file_topology,status)
   call getarg(3,file_output,status)

   write(*,*) "******** GSM SOLVER ********"

   open(unit=1,file=file_output,iostat=ierr)
   if(ierr/=0) then
      write(*,*) "FILE ERROR!; ERROR CODE (OPEN) =",ierr
      write(1,*) "FILE ERROR!; ERROR CODE (OPEN) =",ierr
      stop
   end if

   call input                     ! Input Data (input.f90)

   write(*,*) '**** FINISHED ****'
   write(1,*) '**** FINISHED ****'

   close(unit=1)

   stop

   !--------------------------------
   contains
      subroutine show_help
         write(*,*) "******** GSM SOLVER (VERSION 1.0) ********"
         write(*,*) "  Copyright (c) 2007 Takuichi Hirano (Tokyo Institute of Technology)."
         write(*,*) "  All rights reserved."
         write(*,*)
         write(*,*) "  E-mail: hira@antenna.ee.titech.ac.jp"
         write(*,*) "  http://www-antenna.ee.titech.ac.jp/~hira/"
         write(*,*)
         write(*,*) "Usage: gsm [input S-matrix file] [input topology file] [output file]"
         write(*,*) "  For more information, please see:"
         write(*,*) "  http://www-antenna.ee.titech.ac.jp/~hira/free_software/gsm_solver/"
      end subroutine
end program

!**********************************************************************
! Input Routine
!**********************************************************************

!----------------------------------------------------------------------
! Read Input Parameters
!----------------------------------------------------------------------
subroutine input
   use mod_consts
   use mod_gsm
   implicit none

   intrinsic :: log10
   integer :: i,j,k,n,gport,in_out
   integer :: b1,p1,b2,p2
   real(8) :: mag,pha
   character(len=256) :: text
   integer :: ierr

   write(*,*) '**** INPUT ****'
   write(1,*) '**** INPUT ****'

   ! ******** Read S matrices ********
   ! Table: tab_blk_lport_2_gport and tab_gport_2_blk_lport are made
   ! Identify No. of Total Ports
   write(*,*) "---- READ S MATRICES ----"
   write(1,*) "---- READ S MATRICES ----"

   ! For Array Allocation
   open(unit=10,file=file_s_matrix,iostat=ierr,readonly)
   if(ierr/=0) call file_err    ! Error Trap
   read(unit=10,fmt=*,err=100,end=110) text
   read(unit=10,fmt=*,err=100,end=110) text,n_block
   write(1,*) "NO. OF TOTAL BLOCKS=",n_block
   allocate(n_port(n_block),stat=ierr); call allocate_err_check(ierr)
   m_port=0    ! Maximum number of ports in each block
   n_gsm_total_port=0
   do k=1,n_block
      read(unit=10,fmt=*,err=100,end=110) text,n_port(k)
      n_gsm_total_port=n_gsm_total_port+n_port(k)

      if(n_port(k)>m_port) m_port=n_port(k)    ! renew m_port

      ! read through
      do i=1,n_port(k)
         do j=1,n_port(k)
            read(unit=10,fmt=*,err=100,end=110) text,mag,pha
         end do
      end do
   end do
   write(1,*) "NO. OF TOTAL PORTS=",n_gsm_total_port
   close(unit=10)

   ! In allocate statement, zero clear is guaranteed
   allocate(tab_gport_2_blk_lport(n_gsm_total_port),stat=ierr); call allocate_err_check(ierr)
   allocate(tab_blk_lport_2_gport(n_block,m_port),stat=ierr); call allocate_err_check(ierr)
   allocate(tab_port_connect_info(n_gsm_total_port),stat=ierr); call allocate_err_check(ierr)
   allocate(tab_ex_val_lis(n_gsm_total_port),stat=ierr); call allocate_err_check(ierr)
   allocate(tab_out_lis(n_gsm_total_port),stat=ierr); call allocate_err_check(ierr)
   allocate(tab_unknown_no(n_gsm_total_port,2),stat=ierr); call allocate_err_check(ierr)
   allocate(s_matrix(n_block,m_port,m_port),stat=ierr); call allocate_err_check(ierr)

   ! Read Scattering Matrices
   open(unit=10,file=file_s_matrix,iostat=ierr,readonly)
   if(ierr/=0) call file_err    ! Error Trap
   read(unit=10,fmt=*,err=100,end=110) text
   read(unit=10,fmt=*,err=100,end=110) text,n_block

   n=0    ! n stands for global port no.
   do k=1,n_block
      write(1,*) "-- BLOCK ",k, " --"

      read(unit=10,fmt=*,err=100,end=110) text,n_port(k)

      do i=1,n_port(k)
         ! Make Table
         n=n+1
         ! (Block #, Local Port #) -> Global Port #
         tab_blk_lport_2_gport(k,i)=n
         ! Global Port # -> (Block #, Local Port #)
         tab_gport_2_blk_lport(n).block=k
         tab_gport_2_blk_lport(n).loc_port=i

         do j=1,n_port(k)
            read(unit=10,fmt=*,err=100,end=110) text,mag,pha
            s_matrix(k,i,j)=(10.0d0**(mag/20.0d0))*cdexp(zj*(pha*(pi/180.0d0)))

            if(cdabs(s_matrix(k,i,j))<eps) s_matrix(k,i,j)=dcmplx(eps)
            write(1,*) "S(",i,",",j,")=",20.0d0*log10(cdabs(s_matrix(k,i,j)))," [dB], ", &
                       datan2d(dimag(s_matrix(k,i,j)),dreal(s_matrix(k,i,j)))," [deg]"
         end do
      end do
   end do
   write(1,*) "---- END OF READ S MATRICES ----"
   close(unit=10)

   ! ******** Read Topology ********
   ! Table: tab_port_connect_info and tab_ex_val_lis are made
   write(*,*) "---- READ TOPOLOGY ----"
   write(1,*) "---- READ TOPOLOGY ----"

   ! Zero Clear
   do i=1,n_gsm_total_port
      tab_port_connect_info(i)=0
   end do

   open(unit=10,file=file_topology,iostat=ierr,readonly)
   if(ierr/=0) call file_err    ! Error Trap
   read(unit=10,fmt=*,err=100,end=110) text
   n_output_port=0
   do while(not(eof(10)))
      read(unit=10,fmt='(a2)',advance='no',err=100) text
      select case (text)
         !-------- CM (Comment) --------
         case ("CM")
            read(unit=10,fmt=*,err=100,end=110) text
            write(1,*) "CM"
         !-------- CN (Connect) --------
         case ("CN")
            read(unit=10,fmt=*,err=100) b1,p1,b2,p2
            if((b1>n_block).or.(b2>n_block)) call topology_err(1)          ! Error Trap
            if((p1>n_port(b1)).or.(p2>n_port(b2))) call topology_err(1)    ! Error Trap

            if((tab_port_connect_info(tab_blk_lport_2_gport(b1,p1))/=0).or.&
               (tab_port_connect_info(tab_blk_lport_2_gport(b2,p2))/=0)) call topology_err(2)    ! Error Trap

            tab_port_connect_info(tab_blk_lport_2_gport(b1,p1))=tab_blk_lport_2_gport(b2,p2)
			tab_port_connect_info(tab_blk_lport_2_gport(b2,p2))=tab_blk_lport_2_gport(b1,p1)
            write(1,*) "CN",b1,p1,b2,p2
         !-------- EX (Excitation) --------
         case ("EX")
            read(unit=10,fmt=*,err=100,end=110) b1,p1,mag,pha
            if(b1>n_block) call topology_err(1)       ! Error Trap
            if(p1>n_port(b1)) call topology_err(1)    ! Error Trap

            if((tab_port_connect_info(tab_blk_lport_2_gport(b1,p1))/=0)) call topology_err(2)    ! Error Trap

            tab_port_connect_info(tab_blk_lport_2_gport(b1,p1))=-1
            tab_ex_val_lis(tab_blk_lport_2_gport(b1,p1))=(10.0d0**(mag/20.0d0))*cdexp(zj*(pha*(pi/180.0d0)))
            write(1,*) "EX",b1,p1,mag,pha
         !-------- LD (Matched Load) --------
         case ("LD")
            read(unit=10,fmt=*,err=100,end=110) b1,p1
            if(b1>n_block) call topology_err(1)       ! Error Trap
            if(p1>n_port(b1)) call topology_err(1)    ! Error Trap

            if((tab_port_connect_info(tab_blk_lport_2_gport(b1,p1))/=0)) call topology_err(2)    ! Error Trap

            tab_port_connect_info(tab_blk_lport_2_gport(b1,p1))=0
            write(1,*) "LD",b1,p1
         !-------- OP (Output) --------
         case ("OP")
            read(unit=10,fmt=*,err=100,end=110) b1,p1,in_out
            if(b1>n_block) call topology_err(1)       ! Error Trap
            if(p1>n_port(b1)) call topology_err(1)    ! Error Trap
            n_output_port=n_output_port+1
            tab_out_lis(n_output_port).gport=tab_blk_lport_2_gport(b1,p1)
            tab_out_lis(n_output_port).in_out=in_out
            write(1,*) "OP",b1,p1,in_out
         !-------- ED (END) --------
         case ("ED")
            write(1,*) "ED"
            exit
         !-------- Others --------
         case default
            write(*,*) "UNEXPECTED INPUT (TOPOLOGY FILE)!"
            stop
      end select
   end do
   write(1,*) "---- END OF READ TOPOLOGY ----"
   close(unit=10)

   ! ******** Make Table ********
   ! Table: tab_unknown_no
   n_gsm_matrix_unknown=0
   do k=1,n_block
      do i=1,n_port(k)
         gport=tab_blk_lport_2_gport(k,i)
         select case (tab_port_connect_info(gport))
            !-------- Matched Load --------
            case (0)
               ! input
               tab_unknown_no(gport,1)=0
               ! output
               n_gsm_matrix_unknown=n_gsm_matrix_unknown+1
               tab_unknown_no(gport,2)=n_gsm_matrix_unknown
            !-------- Excite --------
            case (-1)
               ! input
               tab_unknown_no(gport,1)=-1
               ! output
               n_gsm_matrix_unknown=n_gsm_matrix_unknown+1
               tab_unknown_no(gport,2)=n_gsm_matrix_unknown
            !-------- Others --------
            case default
               ! input
               n_gsm_matrix_unknown=n_gsm_matrix_unknown+1
               tab_unknown_no(gport,1)=n_gsm_matrix_unknown
               ! output
               n_gsm_matrix_unknown=n_gsm_matrix_unknown+1
               tab_unknown_no(gport,2)=n_gsm_matrix_unknown
         end select
      end do
   end do
   write(1,*) "NO. OF GSM UNKNOWNS=",n_gsm_matrix_unknown

   call solve    ! Solve GSM System Matrix (solve.f90)

   ! Deallocate Arrays
   deallocate(s_matrix,stat=ierr); call deallocate_err_check(ierr)
   deallocate(tab_unknown_no,stat=ierr); call deallocate_err_check(ierr)
   deallocate(tab_out_lis,stat=ierr); call deallocate_err_check(ierr)
   deallocate(tab_ex_val_lis,stat=ierr); call deallocate_err_check(ierr)
   deallocate(tab_port_connect_info,stat=ierr); call deallocate_err_check(ierr)
   deallocate(tab_blk_lport_2_gport,stat=ierr); call deallocate_err_check(ierr)
   deallocate(tab_gport_2_blk_lport,stat=ierr); call deallocate_err_check(ierr)
   deallocate(n_port,stat=ierr); call deallocate_err_check(ierr)

   return

   !---------------- ERROR TRAP FOR READ ----------------
100 &
   write(*,*) "FILE READ ERROR! (FORMAT ERROR)"
   write(1,*) "FILE READ ERROR! (FORMAT ERROR)"
   stop
110 &
   write(*,*) "FILE READ ERROR! (REACHED EOF UNEXPECTEDLY)"
   write(1,*) "FILE READ ERROR! (REACHED EOF UNEXPECTEDLY)"
   stop

   !--------------------------------
   contains
      !--------
      subroutine file_err
         write(*,*) "FILE ERROR!; ERROR CODE (OPEN) =",ierr
         write(1,*) "FILE ERROR!; ERROR CODE (OPEN) =",ierr
         stop
      end subroutine
      !--------
      subroutine topology_err(ierr)
         integer :: ierr

         write(*,*) "TOPOLOGY DESCRIPTION ERROR!"
         write(1,*) "TOPOLOGY DESCRIPTION ERROR!"

         select case (ierr)
            !-------- OUT OF RANGE --------
            case (1)
               write(*,*) "OUT OF RANGE: BLOCK OR PORT NO."
               write(1,*) "OUT OF RANGE: BLOCK OR PORT NO."
            !-------- DUPLICATED CONNECTION --------
            case (2)
               write(*,*) "CONNECTION OR TERMINATION OVERLAPPING"
               write(1,*) "CONNECTION OR TERMINATION OVERLAPPING"
            !-------- Others --------
            case default
         end select

         stop
      end subroutine
end subroutine

!**********************************************************************
! Solve Routine
!**********************************************************************

!----------------------------------------------------------------------
! Solve GSM Matrix
!----------------------------------------------------------------------
subroutine solve
   use mod_consts
   use mod_gsm
   implicit none

   integer :: i,j,k,idx_i,idx_j,block_i,s_i,s_j
   complex(8) :: c
   integer :: ierr

   allocate(gsm_matrix(n_gsm_matrix_unknown,n_gsm_matrix_unknown),stat=ierr); call allocate_err_check(ierr)
   allocate(gsm_rhs_vec(n_gsm_matrix_unknown),stat=ierr); call allocate_err_check(ierr)

   ! ******** Build GSM Matrix ********

   write(*,*) '**** BUILD GSM MATRIX ****'
   write(1,*) '**** BUILD GSM MATRIX ****'

   ! Zero Clear
   call matrix_zero_clear(gsm_matrix,n_gsm_matrix_unknown,n_gsm_matrix_unknown,n_gsm_matrix_unknown,n_gsm_matrix_unknown)
   call vector_zero_clear(gsm_rhs_vec,n_gsm_matrix_unknown,n_gsm_matrix_unknown)

   do i=1,n_gsm_total_port
      ! ---- Equation for Input ----
      j=tab_port_connect_info(i)      ! Connected Global Port No.
      if(j > 0) then
         idx_i=tab_unknown_no(i,1)    ! Unknown No. of Self Port #i
         idx_j=tab_unknown_no(j,2)    ! Unknown No. of Connected Port #j

         !    Eq. for ##i: c_{i}=c{j}
         ! -> Eq. for ##i: c_{i}-c{j}=0
         gsm_matrix(idx_i,idx_i)=gsm_matrix(idx_i,idx_i)+1.0d0
         gsm_matrix(idx_i,idx_j)=gsm_matrix(idx_i,idx_j)-1.0d0
      end if

      ! ---- Equation for Output ----
      idx_i=tab_unknown_no(i,2)       ! Unknown No. of Self Port #i

      block_i=tab_gport_2_blk_lport(i).block     ! Block No. of Port #i
      s_i=tab_gport_2_blk_lport(i).loc_port      ! Local Port No. of Port #i

      !    for ##i: c_{i}=Sum c{j} S_{ij}
      ! -> for ##i: c_{i}-Sum c{j} S_{ij}=0
      gsm_matrix(idx_i,idx_i)=gsm_matrix(idx_i,idx_i)+1.0d0

      do k=1,n_port(block_i)
         j=tab_blk_lport_2_gport(block_i,k)      ! Global Port No. #k
         idx_j=tab_unknown_no(j,1)               ! Unknown No. of Connected Port #j

         s_j=k

         select case (idx_j)
            case (0)
               ! No Input from Port #j
            case (-1)
               ! Excitation from Port #j
               gsm_rhs_vec(idx_i)=gsm_rhs_vec(idx_i)+tab_ex_val_lis(j)*s_matrix(block_i,s_i,s_j)
            case default
               ! Add to GSM System Matrix
               gsm_matrix(idx_i,idx_j)=gsm_matrix(idx_i,idx_j)-s_matrix(block_i,s_i,s_j)
         end select
      end do
   end do

   ! ******** Solve GSM System Matrix ********
   call gsm_solve_matrix_equation

   ! ******** Output Solution ********
   write(*,*) '**** OUTPUT ****'
   write(1,*) '**** OUTPUT ****'
 
   write(1,*) "(BLOCK, PORT, IN[1] OR OUT[2])=MAGNITUDE [dB], PHASE [deg]"

   do i=1,n_output_port
      select case (tab_unknown_no(tab_out_lis(i).gport,tab_out_lis(i).in_out))
         case (0)
            ! Matched Load (No Input)
            c=dcmplx(eps)
         case (-1)
            ! Excitation (Known)
            c=tab_ex_val_lis(tab_out_lis(i).gport)
         case default
            ! Normal Unknown
            c=gsm_rhs_vec(tab_unknown_no(tab_out_lis(i).gport,tab_out_lis(i).in_out))
      end select

      write(1,*) "(",tab_gport_2_blk_lport(tab_out_lis(i).gport).block,",", &
                 tab_gport_2_blk_lport(tab_out_lis(i).gport).loc_port,",", &
                 tab_out_lis(i).in_out,")="

      if(cdabs(c)<eps) c=dcmplx(eps)
      write(1,*) 20.0d0*log10(cdabs(c))," [dB], ",datan2d(dimag(c),dreal(c))," [deg]"
   end do

   deallocate(gsm_rhs_vec,stat=ierr); call deallocate_err_check(ierr)
   deallocate(gsm_matrix,stat=ierr); call deallocate_err_check(ierr)

   return
end subroutine

!----------------------------------------------------------------------
! Solve GSM System Matrix
!----------------------------------------------------------------------
subroutine gsm_solve_matrix_equation
   use mod_consts
   use mod_gsm
   implicit none

   write(*,*) '**** SOLVE ****'
   write(1,*) '**** SOLVE ****'

   ! **** Solve Matrix Equation ****
   ! A x=b
   ! A -> n_gsm_matrix_unknowns, x -> gsm_rhs_vec, b -> gsm_rhs_vec
   ! --------
   ! Pseudo Inverse
   call gsm_matrix_pseudo_inversion

   return
end subroutine

!----------------------------------------------------------------------
! Pseudo Inverse (Generalized Inverse, Moore-Penrose Inverse)
! Solve matrix equation
! A x=b
!
! Even if det(A)=0,
! the solution x', in a sence of min ||A x'- b||, can be obtained.
!
! Singular value decomposition (SVD) routine of LAPACK is used
! because pseudo inverse routine is not available.
!
! [References]
! G. Strang, Linear Algebra and its Applications, 3rd ed., Thomson Learning, 1988.
!----------------------------------------------------------------------
subroutine gsm_matrix_pseudo_inversion
   use mod_consts
   use mod_gsm
   implicit none

   integer :: i,j

   integer :: info,ierr
   real(8) :: dsum,cond
   real(8),allocatable :: sing_val(:),dwork(:)
   complex(8),allocatable :: u(:,:),vt(:,:),zwork(:)
   complex(8),allocatable :: a_wk(:,:),b_wk1(:),b_wk2(:)
   complex(8) :: zsum

   allocate(sing_val(n_gsm_matrix_unknown),stat=ierr); call allocate_err_check(ierr)
   allocate(dwork(5*n_gsm_matrix_unknown),stat=ierr); call allocate_err_check(ierr)
   allocate(u(n_gsm_matrix_unknown,n_gsm_matrix_unknown),stat=ierr); call allocate_err_check(ierr)
   allocate(vt(n_gsm_matrix_unknown,n_gsm_matrix_unknown),stat=ierr); call allocate_err_check(ierr)
   allocate(zwork(3*n_gsm_matrix_unknown),stat=ierr); call allocate_err_check(ierr)
   allocate(a_wk(n_gsm_matrix_unknown,n_gsm_matrix_unknown),stat=ierr); call allocate_err_check(ierr)
   allocate(b_wk1(n_gsm_matrix_unknown),stat=ierr); call allocate_err_check(ierr)
   allocate(b_wk2(n_gsm_matrix_unknown),stat=ierr); call allocate_err_check(ierr)

   ! Copy because gsm_matrix will be destroied in zgesvd
   do i=1,n_gsm_matrix_unknown
      do j=1,n_gsm_matrix_unknown
         a_wk(i,j)=gsm_matrix(i,j)
      end do
   end do

   ! Singular Value Decomposition (SVD)
   !
   ! Netlib (http://www.netlib.org/)
   ! LAPACK (Linear Algebra PACKage) (http://www.netlib.org/lapack/)
   ! lapack/complex16 (http://www.netlib.org/lapack/complex16/)
   ! "zgesvd.f plus dependencies"
   !
   call zgesvd('A', 'A', n_gsm_matrix_unknown, n_gsm_matrix_unknown, a_wk, n_gsm_matrix_unknown,&
               sing_val, u, n_gsm_matrix_unknown, vt, n_gsm_matrix_unknown, zwork, 3*n_gsm_matrix_unknown, dwork, info)

   if(info.ne.0) then
      write(*,*) "Sub(pseudo_inverse): argument is illegal for ''zgesvd''."
      write(*,*) "error code",info
      stop
   end if

   ! -------- Singular Value --------
   !do i=1,min(m,n)
   !   write(*,*) sing_val(i)
   !end do
   !stop

   ! -------- Condition Number --------
   ! (Condition Number)=(Max Singular Value)/(Min Singular Value)
   cond=sing_val(1)
   do i=1,n_gsm_matrix_unknown
      if(sing_val(i) > eps2) then
         cond=sing_val(i)
      end if
   end do
   cond=sing_val(1)/cond
   write(1,*) "CONDITION NUMBER=",cond

   ! SVD: A    =U    *S    *VT
   !      [m,n] [m,m] [m,n] [n,n]
   ! pseudo_inverse(A)=VT^H*S^(-1)*U^H
   ! M^H means Hermitian conjugate of matrix M 

   ! A x = b
   ! x = pseudo_inverse(A) b
   !   = VT^H*S^(-1)*U^H b

   ! b_wk1:= b
   ! [m,1]   [m,1]
   do i=1,n_gsm_matrix_unknown
      b_wk1(i)=gsm_rhs_vec(i)
   end do

   ! b_wk := U^H   b
   ! [m,1]   [m,m] [m,1]
   do i=1,n_gsm_matrix_unknown
      zsum=dcmplx(0.0d0,0.0d0)
      do j=1,n_gsm_matrix_unknown
         zsum=zsum+dconjg(u(j,i))*b_wk1(j)
      end do
      b_wk2(i)=zsum
   end do

   ! b_wk1:= b_wk2
   ! [m,1]   [m,1]
   do i=1,n_gsm_matrix_unknown
      b_wk1(i)=b_wk2(i)
   end do

   ! b_wk := S^(-1)*U^H b
   ! [n,1]   [n,m]  [m,1]
   do i=1,n_gsm_matrix_unknown
      if(sing_val(i) > eps2) then
         b_wk2(i)=b_wk1(i)/sing_val(i)
      end if
   end do

   ! b_wk1:= b_wk2
   ! [n,1]   [n,1]
   do i=1,n_gsm_matrix_unknown
      b_wk1(i)=b_wk2(i)
   end do

   ! b_wk := VT^H *S^(-1)*U^H b
   ! [n,1]   [n,n] [n,1]
   do i=1,n_gsm_matrix_unknown
      zsum=dcmplx(0.0d0,0.0d0)
      do j=1,n_gsm_matrix_unknown
         zsum=zsum+dconjg(vt(j,i))*b_wk1(j)
      end do
      b_wk2(i)=zsum
   end do

   ! b_wk1:= b_wk2
   ! [n,1]   [n,1]
   do i=1,n_gsm_matrix_unknown
      b_wk1(i)=b_wk2(i)
      gsm_rhs_vec(i)=b_wk2(i)        ! Over write solution to R.H.S. vector
   end do

   ! -------- Calculate Square Error --------
   ! Calculate A x'
   !do i=1,m
   !   zsum=dcmplx(0.0d0)
   !   do j=1,n
   !      zsum=zsum+a(i,j)*b_wk1(j)
   !   end do
   !   b_wk2(i)=zsum
   !end do
   ! Calculate || A x'-b ||
   !dsum=0.0d0
   !do i=1,n
   !   dsum=dsum+cdabs(b_wk2(i)-b(i))**2
   !end do
   !write(*,*) dsum
   !stop

   deallocate(b_wk2,stat=ierr); call deallocate_err_check(ierr)
   deallocate(b_wk1,stat=ierr); call deallocate_err_check(ierr)
   deallocate(a_wk,stat=ierr); call deallocate_err_check(ierr)
   deallocate(zwork,stat=ierr); call deallocate_err_check(ierr)
   deallocate(vt,stat=ierr); call deallocate_err_check(ierr)
   deallocate(u,stat=ierr); call deallocate_err_check(ierr)
   deallocate(dwork,stat=ierr); call deallocate_err_check(ierr)
   deallocate(sing_val,stat=ierr); call deallocate_err_check(ierr)

   return
end subroutine

!----------------------------------------------------------------------
! Zero clear matrix
!
! [Input]
! a: Matrix (Array)
! mm: declared size of rows
! nn: declared size of columns
! m: number of rows
! n: number of columns
!----------------------------------------------------------------------
subroutine matrix_zero_clear(a,mm,nn,m,n)
   implicit none

   integer :: mm,nn,m,n
   complex(8) :: a(mm,nn)

   integer :: i,j

   ! ---- Zero clear matrix ----
   do i=1,m
      do j=1,n
         a(i,j)=dcmplx(0.0d0,0.0d0)
      end do
   end do

   return
end subroutine

!----------------------------------------------------------------------
! Zero clear vector
!
! [Input]
! a: Matrix (Array)
! mm: declared dimension
! m: dimension
!----------------------------------------------------------------------
subroutine vector_zero_clear(b,mm,m)
   implicit none

   integer :: mm,m
   complex(8) :: b(mm)

   integer :: i

   ! ---- Zero clear vector ----
   do i=1,m
      b(i)=dcmplx(0.0d0,0.0d0)
   end do

   return
end subroutine

!----------------------------------------------------------------------
! Allocate error check
!
! [Input]
! ierr: status
!----------------------------------------------------------------------
subroutine allocate_err_check(ierr)
   implicit none
   integer :: ierr

   if(ierr==0) then
      ! successful
      return
   else
      ! fail
      write(*,*) "MEMORY ALLOCATE ERROR!; ERROR CODE (ALLOCATE) =",ierr
      write(1,*) "MEMORY ALLOCATE ERROR!; ERROR CODE (ALLOCATE) =",ierr
      stop
   end if

   return
end subroutine

!----------------------------------------------------------------------
! Deallocate error check
!
! [Input]
! ierr: status
!----------------------------------------------------------------------
subroutine deallocate_err_check(ierr)
   implicit none
   integer :: ierr

   if(ierr==0) then
      ! successful
      return
   else
      ! fail
      write(*,*) "MEMORY DEALLOCATE ERROR!; ERROR CODE (DEALLOCATE) =",ierr
      write(1,*) "MEMORY DEALLOCATE ERROR!; ERROR CODE (DEALLOCATE) =",ierr
      stop
   end if

   return
end subroutine

!
! End of File
!