!********************************************************
!* Solving a complex linear system AX = B by LU decom-  *
!* position with dynamic allocations                    *
!*                                                      *
!*               F90 Release 1.0 By J-P Moreau, Paris   *
!*                        (www.jpmoreau.fr)             *
!* ---------------------------------------------------- *
!* SAMPLE RUN:                                          *
!*                                                      *
!* Input file (test_clu.dat):                           *
!*                                                      *
!* 15                                                   *
!* 0. 0. 0. 0. 6. 0. 1. 0. 1. 0. 4. 0. 2. 0. 1. 0. 2. 0 *. 
!* 2. 0. 0. 0. 3. 0. 0. 0. 5. 0. 0. 1. 1. 1.            *
!* 6. 0. 2. 0. 5. 0. 2. 0. 4. 0. 5. 0. 5. 0. 2. 0. 1. 0.* 
!* 2. 0. 3. 0. 1. 0. 5. 0. 1. 0. 3. 2. 2. 2.            *
!* 6. 0. 2. 0. 5. 0. 6. 0. 3. 0. 6. 0. 5. 0. 0. 0. 0. 0 * 
!* 1. 0. 3. 0. 0. 0. 4. 0. 0. 0. 5. 3. 3. 3.            *
!* 5. 0. 4. 0. 3. 0. 1. 0. 4. 0. 4. 0. 6. 0. 4. 0. 6. 0 *. 
!* 1. 0. 4. 0. 2. 0. 0. 0. 5. 0. 3. 4. 4. 4.            *
!* 3. 0. 4. 0. 6. 0. 4. 0. 2. 0. 0. 0. 4. 0. 6. 0. 5. 0.* 
!* 5. 0. 1. 0. 1. 0. 6. 0. 5. 0. 4. 5. 5. 5.            *
!* 4. 0. 2. 0. 0. 0. 4. 0. 1. 0. 6. 0. 1. 0. 1. 0. 3. 0.* 
!* 6. 0. 1. 0. 1. 0. 1. 0. 1. 0. 4. 6. 6. 6.            *
!* 1. 0. 1. 0. 4. 0. 5. 0. 5. 0. 0. 0. 0. 0. 0. 0. 0. 0.* 
!* 5. 0. 5. 0. 2. 0. 4. 0. 3. 0. 5. 7. 7. 7.            * 
!* 2. 0. 5. 0. 5. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 4. 0.* 
!* 1. 0. 3. 0. 1. 0. 2. 0. 5. 0. 2. 8. 8. 8.            *
!* 4. 0. 4. 0. 6. 0. 2. 0. 3. 0. 2. 0. 0. 0. 5. 0. 5. 0.* 
!* 3  0. 5. 0. 2. 0. 3. 0. 4. 0. 6. 9. 9. 9.            *
!* 3. 0. 2. 0. 0. 0. 0. 0. 3. 0. 4. 0. 0. 0. 6. 0. 5. 0.* 
!* 4. 0. 6. 0. 1. 0. 4. 0. 5. 0. 3.10.10.10.            * 
!* 1. 0. 1. 0. 6. 0. 2. 0. 0. 0. 6. 0. 0. 0. 2. 0. 0. 0.* 
!* 0. 0. 6. 0. 3. 0. 4. 0. 0. 0. 4.11.11.11.            *
!* 6. 0. 3. 0. 2. 0. 3. 0. 6. 0. 2. 0. 0. 0. 2. 0. 3. 0.* 
!* 0. 0. 5. 0. 5. 0. 1. 0. 6. 0. 3.12.12.12.            * 
!* 2. 0. 1. 0. 5. 0. 4. 0. 0. 0. 4. 0. 0. 0. 6. 0. 3. 0.* 
!* 2. 0. 3. 0. 1. 0. 1. 0. 4. 0. 6.13.13.13.            *
!* 3. 0. 0. 0. 5. 0. 1. 0. 1. 0. 6. 0. 0. 0. 4. 0. 3. 0.* 
!* 6. 0. 1. 0. 0. 0. 1. 0. 0. 0. 2.14.14.14.            *
!* 0. 0. 3. 0. 3. 0. 3. 0. 4. 0. 5. 0. 0. 0. 3. 0. 1. 0.* 
!* 2. 0. 4. 0. 6. 0. 4. 0. 1. 0. 4.15.15.15.            *
!*                                                      *
!* Output file (test_clu.lst):                          *
!*                                                      *
!* System solution:                                     *
!*                                                      *
!*  X1 =     -0.057831112295389     0.092401653528214   *
!*  X2 =     -0.885622978210449     1.415035247802734   *
!*  X3 =     -0.106641165912151     0.170389726758003   *
!*  X4 =     -0.421776235103607     0.673907876014709   *
!*  X5 =      0.063331373035908    -0.101189911365509   *
!*  X6 =     -0.098265357315540     0.157006934285164   *
!*  X7 =      0.217855080962181    -0.348085731267929   *
!*  X8 =     -0.698290705680847     1.115718603134155   *
!*  X9 =      0.869955599308014    -1.390002012252808   *
!*  X10=     -0.338608175516129     0.541023075580597   *
!*  X11=     -0.463158965110779     0.740028560161591   *
!*  X12=      0.118877366185188    -0.189940482378006   *
!*  X13=      0.514975011348724    -0.822819411754608   *
!*  X14=      0.013270998373628    -0.021204300224781   *
!*  X15=      0.731169581413269    -1.168251872062683   *
!* ---------------------------------------------------- *
!* Uses: modules BASIS, CLU1                            *
!*                                                      *
!********************************************************
  Program Test_LU

  USE BASIS
  USE CLU1

  COMPLEX, pointer ::  A(:,:)   !complex matrix (n x n)
  COMPLEX, pointer ::  B(:)     !complex vector (n)
  REAL,    pointer ::  temp(:)  !complex temporary vector (n+1)
  integer, pointer ::  INDX(:)  !integer vector (n)

  integer d, rc
  character*12 input, output
  character*8 s
  
! input data file name from screen
  print *,' '
  write(*,50,advance='no'); read *, s

  J=0
  DO I=1,LEN(s)
    IF(s(I:I)<>' ') J=J+1  !J=real length of string s
  ENDDO
      
  input=s(1:J)//'.dat'
  output =s(1:J)//'.lst'

! open input and output files
  open(unit=1,file=input,status='old')
  open(unit=2,file=output,status='unknown')

  read(1,*) n     !size of given linear system

  n2 = 2*n + 2

  !dynamic allocations
  allocate(A(n,n),stat=ialloc)
  allocate(B(n),stat=ialloc)
  allocate(temp(n2),stat=ialloc)
  allocate(INDX(n),stat=ialloc)

!Read data in file #1 and copy to file #2
  write(2,*) ' '
  call WriteHead(2,'  COMPLEX LINEAR SYSTEM TO BE SOLVED:')
  write(2,60) n

  do i=1, n 
    !read one line
    call ReadCVec(1,n2,temp)
	j1=1
    do j=1, n
      A(i,j) = CMPLX(temp(j1),temp(j1+1))
	  j1=j1+2
    end do
    B(i) = CMPLX(temp(n2-1),temp(n2))
    !print one line
    call WriteCVec(2,n2,temp)
  end do
  close(1)


!call LU decomposition routine
  call LUDCMP(A,n,INDX,D,rc)

!call appropriate solver if previous return code is ok
  if (rc.eq.0) then
    call LUBKSB(A,n,INDX,B)
  endif

!print results or error message
  if (rc.eq.1) then
    write(2,*) ' The system matrix is singular, no solution !'
  else
    write(2,*) ' '
    write(2,*) '  System solution:'
    write(2,*)
    if(n<10) then
      do i=1, n 
        write(2,200) i,B(i) 
      end do
    else
      do i=1, 9 
        write(2,200) i,B(i) 
      end do
      do i=10, n 
        write(2,201) i,B(i) 
      end do
    end if
  end if

! section fin
  call WriteEnd(2)
  Close(2)
  Write(*,*) ' Results in file ',output,'.'
  stop
50  format(' Input data file name (without .dat): ')
60  format(/' N = ',I2/)
200 format('    X',I1,' = ',2F22.15)
201 format('    X',I2,'= ',2F22.15)
  stop
  END

  Subroutine ReadCVec (u, n, x) 
  INTEGER u
  REAL x(1:n)
  read(u,*) (x(i),i=1,n)
  return
  end

  Subroutine WriteCVec (u, n, x) 
  INTEGER u
  REAL x(1:n)
  write(u,10) (x(i),i=1,n)
  return
10 format(10F5.1)
  end

! end of file test_clu.f90