home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libblas / ludec / linpack.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  4.1 KB  |  143 lines
  1. /* *****************************************************************************
  2. *
  3. * Copyright 1991, 1992, 1993, 1994, Silicon Graphics, Inc.
  4. * All Rights Reserved.
  5. *
  6. * This is UNPUBLISHED PROPRIETARY SOURCE CODE of Silicon Graphics, Inc.;
  7. * the contents of this file may not be disclosed to third parties, copied or
  8. * duplicated in any form, in whole or in part, without the prior written
  9. * permission of Silicon Graphics, Inc.
  10. *
  11. * RESTRICTED RIGHTS LEGEND:
  12. * Use, duplication or disclosure by the Government is subject to restrictions
  13. * as set forth in subdivision (c)(1)(ii) of the Rights in Technical Data
  14. * and Computer Software clause at DFARS 252.227-7013, and/or in similar or
  15. * successor clauses in the FAR, DOD or NASA FAR Supplement. Unpublished -
  16. * rights reserved under the Copyright Laws of the United States.
  17. *
  18. ***************************************************************************** */
  19.       subroutine matgen(a,lda,n)
  20.       double precision a(lda,1)
  21. c
  22.       init = 1325
  23.       norma = 0.0
  24. c$doacross
  25.       do 35 i = 1,n
  26.           a(i,1) = 0.0
  27.    35 continue
  28.       do 30 j = 1,n
  29.          do 20 i = 1,n
  30.             init = mod(3125*init,65536)
  31.             a(i,j) = (init - 32768.0)/16384.0
  32.    20    continue
  33.    30    continue
  34.       return
  35.       end
  36.       
  37.       subroutine dgefa(a,lda,n,ipvt,info)
  38.       integer lda,n,ipvt(1),info
  39.       double precision a(lda,1)
  40. c
  41. c     dgefa factors a double precision matrix by gaussian elimination.
  42. c
  43. c     dgefa is usually called by dgeco, but it can be called
  44. c     directly with a saving in time if  rcond  is not needed.
  45. c     (time for dgeco) = (1 + 9/n)*(time for dgefa) .
  46. c
  47. c     on entry
  48. c
  49. c        a       double precision(lda, n)
  50. c                the matrix to be factored.
  51. c
  52. c        lda     integer
  53. c                the leading dimension of the array  a .
  54. c
  55. c        n       integer
  56. c                the order of the matrix  a .
  57. c
  58. c     on return
  59. c
  60. c        a       an upper triangular matrix and the multipliers
  61. c                which were used to obtain it.
  62. c                the factorization can be written  a = l*u  where
  63. c                l  is a product of permutation and unit lower
  64. c                triangular matrices and  u  is upper triangular.
  65. c
  66. c        ipvt    integer(n)
  67. c                an integer vector of pivot indices.
  68. c
  69. c        info    integer
  70. c                = 0  normal value.
  71. c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
  72. c                     condition for this subroutine, but it does
  73. c                     indicate that dgesl or dgedi will divide by zero
  74. c                     if called.  use  rcond  in dgeco for a reliable
  75. c                     indication of singularity.
  76. c
  77. c     linpack. this version dated 08/14/78 .
  78. c     cleve moler, university of new mexico, argonne national lab.
  79. c
  80. c     subroutines and functions
  81. c
  82. c     blas daxpy,dscal,idamax
  83. c
  84. c     internal variables
  85. c
  86.       double precision t
  87.       integer idamax,j,k,kp1,l,nm1
  88. c
  89. c
  90. c     gaussian elimination with partial pivoting
  91. c
  92.       info = 0
  93.       nm1 = n - 1
  94.       if (nm1 .lt. 1) go to 70
  95.       do 60 k = 1, nm1
  96.          kp1 = k + 1
  97. c
  98. c        find l = pivot index
  99. c
  100.          l = idamax(n-k+1,a(k,k),1) + k - 1
  101.          ipvt(k) = l
  102. c
  103. c        zero pivot implies this column already triangularized
  104. c
  105.          if (a(l,k) .eq. 0.0d0) go to 40
  106. c
  107. c           interchange if necessary
  108. c
  109.             if (l .eq. k) go to 10
  110.                t = a(l,k)
  111.                a(l,k) = a(k,k)
  112.                a(k,k) = t
  113.    10       continue
  114. c
  115. c           compute multipliers
  116. c
  117.             t = -1.0d0/a(k,k)
  118.             call dscal(n-k,t,a(k+1,k),1)
  119. c
  120. c           row elimination with column indexing
  121. c
  122. c$doacross local(j,t)
  123.             do 30 j = kp1, n
  124.                t = a(l,j)
  125.                if (l .eq. k) go to 20
  126.                   a(l,j) = a(k,j)
  127.                   a(k,j) = t
  128.    20          continue
  129.                call daxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
  130.    30       continue
  131.          go to 50
  132.    40    continue
  133.             info = k
  134.    50    continue
  135.     if( mod(k,32) .eq. 0) call update(k)
  136.    60 continue
  137.    70 continue
  138.       ipvt(n) = n
  139.       if (a(n,n) .eq. 0.0d0) info = n
  140.     call update(n)
  141.       return
  142.       end
  143.