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 / blas / dgemv.f < prev    next >
Text File  |  1996-09-28  |  10KB  |  322 lines
  1. *
  2. ************************************************************************
  3. *
  4. *     File of the DOUBLE PRECISION  Level-2 BLAS.
  5. *     ===========================================
  6. *
  7. *     SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX,
  8. *    $                   BETA, Y, INCY )
  9. *
  10. *     SUBROUTINE DGBMV ( TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX,
  11. *    $                   BETA, Y, INCY )
  12. *
  13. *     SUBROUTINE DSYMV ( UPLO, N, ALPHA, A, LDA, X, INCX,
  14. *    $                   BETA, Y, INCY )
  15. *
  16. *     SUBROUTINE DSBMV ( UPLO, N, K, ALPHA, A, LDA, X, INCX,
  17. *    $                   BETA, Y, INCY )
  18. *
  19. *     SUBROUTINE DSPMV ( UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY )
  20. *
  21. *     SUBROUTINE DTRMV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX )
  22. *
  23. *     SUBROUTINE DTBMV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX )
  24. *
  25. *     SUBROUTINE DTPMV ( UPLO, TRANS, DIAG, N, AP, X, INCX )
  26. *
  27. *     SUBROUTINE DTRSV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX )
  28. *
  29. *     SUBROUTINE DTBSV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX )
  30. *
  31. *     SUBROUTINE DTPSV ( UPLO, TRANS, DIAG, N, AP, X, INCX )
  32. *
  33. *     SUBROUTINE DGER  ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA )
  34. *
  35. *     SUBROUTINE DSYR  ( UPLO, N, ALPHA, X, INCX, A, LDA )
  36. *
  37. *     SUBROUTINE DSPR  ( UPLO, N, ALPHA, X, INCX, AP )
  38. *
  39. *     SUBROUTINE DSYR2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA )
  40. *
  41. *     SUBROUTINE DSPR2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, AP )
  42. *
  43. *     See:
  44. *
  45. *        Dongarra J. J., Du Croz J. J., Hammarling S.  and Hanson R. J..
  46. *        An  extended  set of Fortran  Basic Linear Algebra Subprograms.
  47. *
  48. *        Technical  Memoranda  Nos. 41 (revision 3) and 81,  Mathematics
  49. *        and  Computer Science  Division,  Argonne  National Laboratory,
  50. *        9700 South Cass Avenue, Argonne, Illinois 60439, US.
  51. *
  52. *        Or
  53. *
  54. *        NAG  Technical Reports TR3/87 and TR4/87,  Numerical Algorithms
  55. *        Group  Ltd.,  NAG  Central  Office,  256  Banbury  Road, Oxford
  56. *        OX2 7DE, UK,  and  Numerical Algorithms Group Inc.,  1101  31st
  57. *        Street,  Suite 100,  Downers Grove,  Illinois 60515-1263,  USA.
  58. *
  59. ************************************************************************
  60. *
  61.       SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX,
  62.      $                   BETA, Y, INCY )
  63. *     .. Scalar Arguments ..
  64.       DOUBLE PRECISION   ALPHA, BETA
  65.       INTEGER            INCX, INCY, LDA, M, N
  66.       CHARACTER*1        TRANS
  67. *     .. Array Arguments ..
  68.       DOUBLE PRECISION   A( LDA, * ), X( * ), Y( * )
  69. *     ..
  70. *
  71. *  Purpose
  72. *  =======
  73. *
  74. *  DGEMV  performs one of the matrix-vector operations
  75. *
  76. *     y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,
  77. *
  78. *  where alpha and beta are scalars, x and y are vectors and A is an
  79. *  m by n matrix.
  80. *
  81. *  Parameters
  82. *  ==========
  83. *
  84. *  TRANS  - CHARACTER*1.
  85. *           On entry, TRANS specifies the operation to be performed as
  86. *           follows:
  87. *
  88. *              TRANS = 'N' or 'n'   y := alpha*A*x + beta*y.
  89. *
  90. *              TRANS = 'T' or 't'   y := alpha*A'*x + beta*y.
  91. *
  92. *              TRANS = 'C' or 'c'   y := alpha*A'*x + beta*y.
  93. *
  94. *           Unchanged on exit.
  95. *
  96. *  M      - INTEGER.
  97. *           On entry, M specifies the number of rows of the matrix A.
  98. *           M must be at least zero.
  99. *           Unchanged on exit.
  100. *
  101. *  N      - INTEGER.
  102. *           On entry, N specifies the number of columns of the matrix A.
  103. *           N must be at least zero.
  104. *           Unchanged on exit.
  105. *
  106. *  ALPHA  - DOUBLE PRECISION.
  107. *           On entry, ALPHA specifies the scalar alpha.
  108. *           Unchanged on exit.
  109. *
  110. *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
  111. *           Before entry, the leading m by n part of the array A must
  112. *           contain the matrix of coefficients.
  113. *           Unchanged on exit.
  114. *
  115. *  LDA    - INTEGER.
  116. *           On entry, LDA specifies the first dimension of A as declared
  117. *           in the calling (sub) program. LDA must be at least
  118. *           max( 1, m ).
  119. *           Unchanged on exit.
  120. *
  121. *  X      - DOUBLE PRECISION array of DIMENSION at least
  122. *           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
  123. *           and at least
  124. *           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
  125. *           Before entry, the incremented array X must contain the
  126. *           vector x.
  127. *           Unchanged on exit.
  128. *
  129. *  INCX   - INTEGER.
  130. *           On entry, INCX specifies the increment for the elements of
  131. *           X. INCX must not be zero.
  132. *           Unchanged on exit.
  133. *
  134. *  BETA   - DOUBLE PRECISION.
  135. *           On entry, BETA specifies the scalar beta. When BETA is
  136. *           supplied as zero then Y need not be set on input.
  137. *           Unchanged on exit.
  138. *
  139. *  Y      - DOUBLE PRECISION array of DIMENSION at least
  140. *           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
  141. *           and at least
  142. *           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
  143. *           Before entry with BETA non-zero, the incremented array Y
  144. *           must contain the vector y. On exit, Y is overwritten by the
  145. *           updated vector y.
  146. *
  147. *  INCY   - INTEGER.
  148. *           On entry, INCY specifies the increment for the elements of
  149. *           Y. INCY must not be zero.
  150. *           Unchanged on exit.
  151. *
  152. *
  153. *  Level 2 Blas routine.
  154. *
  155. *  -- Written on 22-October-1986.
  156. *     Jack Dongarra, Argonne National Lab.
  157. *     Jeremy Du Croz, Nag Central Office.
  158. *     Sven Hammarling, Nag Central Office.
  159. *     Richard Hanson, Sandia National Labs.
  160. *
  161. *
  162. *     .. Parameters ..
  163.       DOUBLE PRECISION   ONE         , ZERO
  164.       PARAMETER        ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  165. *     .. Local Scalars ..
  166.       DOUBLE PRECISION   TEMP
  167.       INTEGER            I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY
  168. *     .. External Functions ..
  169.       LOGICAL            LSAME
  170.       EXTERNAL           LSAME
  171. *     .. External Subroutines ..
  172.       EXTERNAL           XERBLA
  173. *     .. Intrinsic Functions ..
  174.       INTRINSIC          MAX
  175. *     ..
  176. *     .. Executable Statements ..
  177. *
  178. *     Test the input parameters.
  179. *
  180.       INFO = 0
  181.       IF     ( .NOT.LSAME( TRANS, 'N' ).AND.
  182.      $         .NOT.LSAME( TRANS, 'T' ).AND.
  183.      $         .NOT.LSAME( TRANS, 'C' )      )THEN
  184.          INFO = 1
  185.       ELSE IF( M.LT.0 )THEN
  186.          INFO = 2
  187.       ELSE IF( N.LT.0 )THEN
  188.          INFO = 3
  189.       ELSE IF( LDA.LT.MAX( 1, M ) )THEN
  190.          INFO = 6
  191.       ELSE IF( INCX.EQ.0 )THEN
  192.          INFO = 8
  193.       ELSE IF( INCY.EQ.0 )THEN
  194.          INFO = 11
  195.       END IF
  196.       IF( INFO.NE.0 )THEN
  197.          CALL XERBLA( 'DGEMV ', INFO )
  198.          RETURN
  199.       END IF
  200. *
  201. *     Quick return if possible.
  202. *
  203.       IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
  204.      $    ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
  205.      $   RETURN
  206. *
  207. *     Set  LENX  and  LENY, the lengths of the vectors x and y, and set
  208. *     up the start points in  X  and  Y.
  209. *
  210.       IF( LSAME( TRANS, 'N' ) )THEN
  211.          LENX = N
  212.          LENY = M
  213.       ELSE
  214.          LENX = M
  215.          LENY = N
  216.       END IF
  217.       IF( INCX.GT.0 )THEN
  218.          KX = 1
  219.       ELSE
  220.          KX = 1 - ( LENX - 1 )*INCX
  221.       END IF
  222.       IF( INCY.GT.0 )THEN
  223.          KY = 1
  224.       ELSE
  225.          KY = 1 - ( LENY - 1 )*INCY
  226.       END IF
  227. *
  228. *     Start the operations. In this version the elements of A are
  229. *     accessed sequentially with one pass through A.
  230. *
  231. *     First form  y := beta*y.
  232. *
  233.       IF( BETA.NE.ONE )THEN
  234.          IF( INCY.EQ.1 )THEN
  235.             IF( BETA.EQ.ZERO )THEN
  236.                DO 10, I = 1, LENY
  237.                   Y( I ) = ZERO
  238.    10          CONTINUE
  239.             ELSE
  240.                DO 20, I = 1, LENY
  241.                   Y( I ) = BETA*Y( I )
  242.    20          CONTINUE
  243.             END IF
  244.          ELSE
  245.             IY = KY
  246.             IF( BETA.EQ.ZERO )THEN
  247.                DO 30, I = 1, LENY
  248.                   Y( IY ) = ZERO
  249.                   IY      = IY   + INCY
  250.    30          CONTINUE
  251.             ELSE
  252.                DO 40, I = 1, LENY
  253.                   Y( IY ) = BETA*Y( IY )
  254.                   IY      = IY           + INCY
  255.    40          CONTINUE
  256.             END IF
  257.          END IF
  258.       END IF
  259.       IF( ALPHA.EQ.ZERO )
  260.      $   RETURN
  261.       IF( LSAME( TRANS, 'N' ) )THEN
  262. *
  263. *        Form  y := alpha*A*x + y.
  264. *
  265.          JX = KX
  266.          IF( INCY.EQ.1 )THEN
  267.             DO 60, J = 1, N
  268.                IF( X( JX ).NE.ZERO )THEN
  269.                   TEMP = ALPHA*X( JX )
  270.                   DO 50, I = 1, M
  271.                      Y( I ) = Y( I ) + TEMP*A( I, J )
  272.    50             CONTINUE
  273.                END IF
  274.                JX = JX + INCX
  275.    60       CONTINUE
  276.          ELSE
  277.             DO 80, J = 1, N
  278.                IF( X( JX ).NE.ZERO )THEN
  279.                   TEMP = ALPHA*X( JX )
  280.                   IY   = KY
  281.                   DO 70, I = 1, M
  282.                      Y( IY ) = Y( IY ) + TEMP*A( I, J )
  283.                      IY      = IY      + INCY
  284.    70             CONTINUE
  285.                END IF
  286.                JX = JX + INCX
  287.    80       CONTINUE
  288.          END IF
  289.       ELSE
  290. *
  291. *        Form  y := alpha*A'*x + y.
  292. *
  293.          JY = KY
  294.          IF( INCX.EQ.1 )THEN
  295.             DO 100, J = 1, N
  296.                TEMP = ZERO
  297.                DO 90, I = 1, M
  298.                   TEMP = TEMP + A( I, J )*X( I )
  299.    90          CONTINUE
  300.                Y( JY ) = Y( JY ) + ALPHA*TEMP
  301.                JY      = JY      + INCY
  302.   100       CONTINUE
  303.          ELSE
  304.             DO 120, J = 1, N
  305.                TEMP = ZERO
  306.                IX   = KX
  307.                DO 110, I = 1, M
  308.                   TEMP = TEMP + A( I, J )*X( IX )
  309.                   IX   = IX   + INCX
  310.   110          CONTINUE
  311.                Y( JY ) = Y( JY ) + ALPHA*TEMP
  312.                JY      = JY      + INCY
  313.   120       CONTINUE
  314.          END IF
  315.       END IF
  316. *
  317.       RETURN
  318. *
  319. *     End of DGEMV .
  320. *
  321.       END
  322.