home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / linpack / zgesl.f < prev   
Text File  |  1996-09-28  |  3KB  |  123 lines
  1.       subroutine zgesl(a,lda,n,ipvt,b,job)
  2.       integer lda,n,ipvt(1),job
  3.       complex*16 a(lda,1),b(1)
  4. c
  5. c     zgesl solves the complex*16 system
  6. c     a * x = b  or  ctrans(a) * x = b
  7. c     using the factors computed by zgeco or zgefa.
  8. c
  9. c     on entry
  10. c
  11. c        a       complex*16(lda, n)
  12. c                the output from zgeco or zgefa.
  13. c
  14. c        lda     integer
  15. c                the leading dimension of the array  a .
  16. c
  17. c        n       integer
  18. c                the order of the matrix  a .
  19. c
  20. c        ipvt    integer(n)
  21. c                the pivot vector from zgeco or zgefa.
  22. c
  23. c        b       complex*16(n)
  24. c                the right hand side vector.
  25. c
  26. c        job     integer
  27. c                = 0         to solve  a*x = b ,
  28. c                = nonzero   to solve  ctrans(a)*x = b  where
  29. c                            ctrans(a)  is the conjugate transpose.
  30. c
  31. c     on return
  32. c
  33. c        b       the solution vector  x .
  34. c
  35. c     error condition
  36. c
  37. c        a division by zero will occur if the input factor contains a
  38. c        zero on the diagonal.  technically this indicates singularity
  39. c        but it is often caused by improper arguments or improper
  40. c        setting of lda .  it will not occur if the subroutines are
  41. c        called correctly and if zgeco has set rcond .gt. 0.0
  42. c        or zgefa has set info .eq. 0 .
  43. c
  44. c     to compute  inverse(a) * c  where  c  is a matrix
  45. c     with  p  columns
  46. c           call zgeco(a,lda,n,ipvt,rcond,z)
  47. c           if (rcond is too small) go to ...
  48. c           do 10 j = 1, p
  49. c              call zgesl(a,lda,n,ipvt,c(1,j),0)
  50. c        10 continue
  51. c
  52. c     linpack. this version dated 08/14/78 .
  53. c     cleve moler, university of new mexico, argonne national lab.
  54. c
  55. c     subroutines and functions
  56. c
  57. c     blas zaxpy,zdotc
  58. c     fortran dconjg
  59. c
  60. c     internal variables
  61. c
  62.       complex*16 zdotc,t
  63.       integer k,kb,l,nm1
  64. c     double precision dreal,dimag
  65. c     complex*16 zdumr,zdumi
  66. c     dreal(zdumr) = zdumr
  67. c     dimag(zdumi) = (0.0d0,-1.0d0)*zdumi
  68. c
  69.       nm1 = n - 1
  70.       if (job .ne. 0) go to 50
  71. c
  72. c        job = 0 , solve  a * x = b
  73. c        first solve  l*y = b
  74. c
  75.          if (nm1 .lt. 1) go to 30
  76.          do 20 k = 1, nm1
  77.             l = ipvt(k)
  78.             t = b(l)
  79.             if (l .eq. k) go to 10
  80.                b(l) = b(k)
  81.                b(k) = t
  82.    10       continue
  83.             call zaxpy(n-k,t,a(k+1,k),1,b(k+1),1)
  84.    20    continue
  85.    30    continue
  86. c
  87. c        now solve  u*x = y
  88. c
  89.          do 40 kb = 1, n
  90.             k = n + 1 - kb
  91.             b(k) = b(k)/a(k,k)
  92.             t = -b(k)
  93.             call zaxpy(k-1,t,a(1,k),1,b(1),1)
  94.    40    continue
  95.       go to 100
  96.    50 continue
  97. c
  98. c        job = nonzero, solve  ctrans(a) * x = b
  99. c        first solve  ctrans(u)*y = b
  100. c
  101.          do 60 k = 1, n
  102.             t = zdotc(k-1,a(1,k),1,b(1),1)
  103.             b(k) = (b(k) - t)/dconjg(a(k,k))
  104.    60    continue
  105. c
  106. c        now solve ctrans(l)*x = y
  107. c
  108.          if (nm1 .lt. 1) go to 90
  109.          do 80 kb = 1, nm1
  110.             k = n - kb
  111.             b(k) = b(k) + zdotc(n-k,a(k+1,k),1,b(k+1),1)
  112.             l = ipvt(k)
  113.             if (l .eq. k) go to 70
  114.                t = b(l)
  115.                b(l) = b(k)
  116.                b(k) = t
  117.    70       continue
  118.    80    continue
  119.    90    continue
  120.   100 continue
  121.       return
  122.       end
  123.