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 / dsyrk.f < prev    next >
Text File  |  1996-09-28  |  9KB  |  298 lines
  1. *
  2. ************************************************************************
  3. *
  4.       SUBROUTINE DSYRK ( UPLO, TRANS, N, K, ALPHA, A, LDA,
  5.      $                   BETA, C, LDC )
  6. *     .. Scalar Arguments ..
  7.       CHARACTER*1        UPLO, TRANS
  8.       INTEGER            N, K, LDA, LDC
  9.       DOUBLE PRECISION   ALPHA, BETA
  10. *     .. Array Arguments ..
  11.       DOUBLE PRECISION   A( LDA, * ), C( LDC, * )
  12. *     ..
  13. *
  14. *  Purpose
  15. *  =======
  16. *
  17. *  DSYRK  performs one of the symmetric rank k operations
  18. *
  19. *     C := alpha*A*A' + beta*C,
  20. *
  21. *  or
  22. *
  23. *     C := alpha*A'*A + beta*C,
  24. *
  25. *  where  alpha and beta  are scalars, C is an  n by n  symmetric matrix
  26. *  and  A  is an  n by k  matrix in the first case and a  k by n  matrix
  27. *  in the second case.
  28. *
  29. *  Parameters
  30. *  ==========
  31. *
  32. *  UPLO   - CHARACTER*1.
  33. *           On  entry,   UPLO  specifies  whether  the  upper  or  lower
  34. *           triangular  part  of the  array  C  is to be  referenced  as
  35. *           follows:
  36. *
  37. *              UPLO = 'U' or 'u'   Only the  upper triangular part of  C
  38. *                                  is to be referenced.
  39. *
  40. *              UPLO = 'L' or 'l'   Only the  lower triangular part of  C
  41. *                                  is to be referenced.
  42. *
  43. *           Unchanged on exit.
  44. *
  45. *  TRANS  - CHARACTER*1.
  46. *           On entry,  TRANS  specifies the operation to be performed as
  47. *           follows:
  48. *
  49. *              TRANS = 'N' or 'n'   C := alpha*A*A' + beta*C.
  50. *
  51. *              TRANS = 'T' or 't'   C := alpha*A'*A + beta*C.
  52. *
  53. *              TRANS = 'C' or 'c'   C := alpha*A'*A + beta*C.
  54. *
  55. *           Unchanged on exit.
  56. *
  57. *  N      - INTEGER.
  58. *           On entry,  N specifies the order of the matrix C.  N must be
  59. *           at least zero.
  60. *           Unchanged on exit.
  61. *
  62. *  K      - INTEGER.
  63. *           On entry with  TRANS = 'N' or 'n',  K  specifies  the number
  64. *           of  columns   of  the   matrix   A,   and  on   entry   with
  65. *           TRANS = 'T' or 't' or 'C' or 'c',  K  specifies  the  number
  66. *           of rows of the matrix  A.  K must be at least zero.
  67. *           Unchanged on exit.
  68. *
  69. *  ALPHA  - DOUBLE PRECISION.
  70. *           On entry, ALPHA specifies the scalar alpha.
  71. *           Unchanged on exit.
  72. *
  73. *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
  74. *           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
  75. *           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
  76. *           part of the array  A  must contain the matrix  A,  otherwise
  77. *           the leading  k by n  part of the array  A  must contain  the
  78. *           matrix A.
  79. *           Unchanged on exit.
  80. *
  81. *  LDA    - INTEGER.
  82. *           On entry, LDA specifies the first dimension of A as declared
  83. *           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
  84. *           then  LDA must be at least  max( 1, n ), otherwise  LDA must
  85. *           be at least  max( 1, k ).
  86. *           Unchanged on exit.
  87. *
  88. *  BETA   - DOUBLE PRECISION.
  89. *           On entry, BETA specifies the scalar beta.
  90. *           Unchanged on exit.
  91. *
  92. *  C      - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
  93. *           Before entry  with  UPLO = 'U' or 'u',  the leading  n by n
  94. *           upper triangular part of the array C must contain the upper
  95. *           triangular part  of the  symmetric matrix  and the strictly
  96. *           lower triangular part of C is not referenced.  On exit, the
  97. *           upper triangular part of the array  C is overwritten by the
  98. *           upper triangular part of the updated matrix.
  99. *           Before entry  with  UPLO = 'L' or 'l',  the leading  n by n
  100. *           lower triangular part of the array C must contain the lower
  101. *           triangular part  of the  symmetric matrix  and the strictly
  102. *           upper triangular part of C is not referenced.  On exit, the
  103. *           lower triangular part of the array  C is overwritten by the
  104. *           lower triangular part of the updated matrix.
  105. *
  106. *  LDC    - INTEGER.
  107. *           On entry, LDC specifies the first dimension of C as declared
  108. *           in  the  calling  (sub)  program.   LDC  must  be  at  least
  109. *           max( 1, n ).
  110. *           Unchanged on exit.
  111. *
  112. *
  113. *  Level 3 Blas routine.
  114. *
  115. *  -- Written on 8-February-1989.
  116. *     Jack Dongarra, Argonne National Laboratory.
  117. *     Iain Duff, AERE Harwell.
  118. *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
  119. *     Sven Hammarling, Numerical Algorithms Group Ltd.
  120. *
  121. *
  122. *     .. External Functions ..
  123.       LOGICAL            LSAME
  124.       EXTERNAL           LSAME
  125. *     .. External Subroutines ..
  126.       EXTERNAL           XERBLA
  127. *     .. Intrinsic Functions ..
  128.       INTRINSIC          MAX
  129. *     .. Local Scalars ..
  130.       LOGICAL            UPPER
  131.       INTEGER            I, INFO, J, L, NROWA
  132.       DOUBLE PRECISION   TEMP
  133. *     .. Parameters ..
  134.       DOUBLE PRECISION   ONE ,         ZERO
  135.       PARAMETER        ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  136. *     ..
  137. *     .. Executable Statements ..
  138. *
  139. *     Test the input parameters.
  140. *
  141.       IF( LSAME( TRANS, 'N' ) )THEN
  142.          NROWA = N
  143.       ELSE
  144.          NROWA = K
  145.       END IF
  146.       UPPER = LSAME( UPLO, 'U' )
  147. *
  148.       INFO = 0
  149.       IF(      ( .NOT.UPPER               ).AND.
  150.      $         ( .NOT.LSAME( UPLO , 'L' ) )      )THEN
  151.          INFO = 1
  152.       ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ).AND.
  153.      $         ( .NOT.LSAME( TRANS, 'T' ) ).AND.
  154.      $         ( .NOT.LSAME( TRANS, 'C' ) )      )THEN
  155.          INFO = 2
  156.       ELSE IF( N  .LT.0               )THEN
  157.          INFO = 3
  158.       ELSE IF( K  .LT.0               )THEN
  159.          INFO = 4
  160.       ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
  161.          INFO = 7
  162.       ELSE IF( LDC.LT.MAX( 1, N     ) )THEN
  163.          INFO = 10
  164.       END IF
  165.       IF( INFO.NE.0 )THEN
  166.          CALL XERBLA( 'DSYRK ', INFO )
  167.          RETURN
  168.       END IF
  169. *
  170. *     Quick return if possible.
  171. *
  172.       IF( ( N.EQ.0 ).OR.
  173.      $    ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
  174.      $   RETURN
  175. *
  176. *     And when  alpha.eq.zero.
  177. *
  178.       IF( ALPHA.EQ.ZERO )THEN
  179.          IF( UPPER )THEN
  180.             IF( BETA.EQ.ZERO )THEN
  181.                DO 20, J = 1, N
  182.                   DO 10, I = 1, J
  183.                      C( I, J ) = ZERO
  184.    10             CONTINUE
  185.    20          CONTINUE
  186.             ELSE
  187.                DO 40, J = 1, N
  188.                   DO 30, I = 1, J
  189.                      C( I, J ) = BETA*C( I, J )
  190.    30             CONTINUE
  191.    40          CONTINUE
  192.             END IF
  193.          ELSE
  194.             IF( BETA.EQ.ZERO )THEN
  195.                DO 60, J = 1, N
  196.                   DO 50, I = J, N
  197.                      C( I, J ) = ZERO
  198.    50             CONTINUE
  199.    60          CONTINUE
  200.             ELSE
  201.                DO 80, J = 1, N
  202.                   DO 70, I = J, N
  203.                      C( I, J ) = BETA*C( I, J )
  204.    70             CONTINUE
  205.    80          CONTINUE
  206.             END IF
  207.          END IF
  208.          RETURN
  209.       END IF
  210. *
  211. *     Start the operations.
  212. *
  213.       IF( LSAME( TRANS, 'N' ) )THEN
  214. *
  215. *        Form  C := alpha*A*A' + beta*C.
  216. *
  217.          IF( UPPER )THEN
  218.             DO 130, J = 1, N
  219.                IF( BETA.EQ.ZERO )THEN
  220.                   DO 90, I = 1, J
  221.                      C( I, J ) = ZERO
  222.    90             CONTINUE
  223.                ELSE IF( BETA.NE.ONE )THEN
  224.                   DO 100, I = 1, J
  225.                      C( I, J ) = BETA*C( I, J )
  226.   100             CONTINUE
  227.                END IF
  228.                DO 120, L = 1, K
  229.                   IF( A( J, L ).NE.ZERO )THEN
  230.                      TEMP = ALPHA*A( J, L )
  231.                      DO 110, I = 1, J
  232.                         C( I, J ) = C( I, J ) + TEMP*A( I, L )
  233.   110                CONTINUE
  234.                   END IF
  235.   120          CONTINUE
  236.   130       CONTINUE
  237.          ELSE
  238.             DO 180, J = 1, N
  239.                IF( BETA.EQ.ZERO )THEN
  240.                   DO 140, I = J, N
  241.                      C( I, J ) = ZERO
  242.   140             CONTINUE
  243.                ELSE IF( BETA.NE.ONE )THEN
  244.                   DO 150, I = J, N
  245.                      C( I, J ) = BETA*C( I, J )
  246.   150             CONTINUE
  247.                END IF
  248.                DO 170, L = 1, K
  249.                   IF( A( J, L ).NE.ZERO )THEN
  250.                      TEMP      = ALPHA*A( J, L )
  251.                      DO 160, I = J, N
  252.                         C( I, J ) = C( I, J ) + TEMP*A( I, L )
  253.   160                CONTINUE
  254.                   END IF
  255.   170          CONTINUE
  256.   180       CONTINUE
  257.          END IF
  258.       ELSE
  259. *
  260. *        Form  C := alpha*A'*A + beta*C.
  261. *
  262.          IF( UPPER )THEN
  263.             DO 210, J = 1, N
  264.                DO 200, I = 1, J
  265.                   TEMP = ZERO
  266.                   DO 190, L = 1, K
  267.                      TEMP = TEMP + A( L, I )*A( L, J )
  268.   190             CONTINUE
  269.                   IF( BETA.EQ.ZERO )THEN
  270.                      C( I, J ) = ALPHA*TEMP
  271.                   ELSE
  272.                      C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
  273.                   END IF
  274.   200          CONTINUE
  275.   210       CONTINUE
  276.          ELSE
  277.             DO 240, J = 1, N
  278.                DO 230, I = J, N
  279.                   TEMP = ZERO
  280.                   DO 220, L = 1, K
  281.                      TEMP = TEMP + A( L, I )*A( L, J )
  282.   220             CONTINUE
  283.                   IF( BETA.EQ.ZERO )THEN
  284.                      C( I, J ) = ALPHA*TEMP
  285.                   ELSE
  286.                      C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
  287.                   END IF
  288.   230          CONTINUE
  289.   240       CONTINUE
  290.          END IF
  291.       END IF
  292. *
  293.       RETURN
  294. *
  295. *     End of DSYRK .
  296. *
  297.       END
  298.