SRC/ssfrk.f

   1 *> \brief \b SSFRK performs a symmetric rank-k operation for matrix in RFP format.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SSFRK + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssfrk.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssfrk.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssfrk.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SSFRK( TRANSR, UPLO, TRANS, N, K, ALPHA, A, LDA, BETA,
  22 *                         C )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       REAL               ALPHA, BETA
  26 *       INTEGER            K, LDA, N
  27 *       CHARACTER          TRANS, TRANSR, UPLO
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       REAL               A( LDA, * ), C( * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> Level 3 BLAS like routine for C in RFP Format.
  40 *>
  41 *> SSFRK performs one of the symmetric rank--k operations
  42 *>
  43 *>    C := alpha*A*A**T + beta*C,
  44 *>
  45 *> or
  46 *>
  47 *>    C := alpha*A**T*A + beta*C,
  48 *>
  49 *> where alpha and beta are real scalars, C is an n--by--n symmetric
  50 *> matrix and A is an n--by--k matrix in the first case and a k--by--n
  51 *> matrix in the second case.
  52 *> \endverbatim
  53 *
  54 *  Arguments:
  55 *  ==========
  56 *
  57 *> \param[in] TRANSR
  58 *> \verbatim
  59 *>          TRANSR is CHARACTER*1
  60 *>          = 'N':  The Normal Form of RFP A is stored;
  61 *>          = 'T':  The Transpose Form of RFP A is stored.
  62 *> \endverbatim
  63 *>
  64 *> \param[in] UPLO
  65 *> \verbatim
  66 *>          UPLO is CHARACTER*1
  67 *>           On  entry, UPLO specifies whether the upper or lower
  68 *>           triangular part of the array C is to be referenced as
  69 *>           follows:
  70 *>
  71 *>              UPLO = 'U' or 'u'   Only the upper triangular part of C
  72 *>                                  is to be referenced.
  73 *>
  74 *>              UPLO = 'L' or 'l'   Only the lower triangular part of C
  75 *>                                  is to be referenced.
  76 *>
  77 *>           Unchanged on exit.
  78 *> \endverbatim
  79 *>
  80 *> \param[in] TRANS
  81 *> \verbatim
  82 *>          TRANS is CHARACTER*1
  83 *>           On entry, TRANS specifies the operation to be performed as
  84 *>           follows:
  85 *>
  86 *>              TRANS = 'N' or 'n'   C := alpha*A*A**T + beta*C.
  87 *>
  88 *>              TRANS = 'T' or 't'   C := alpha*A**T*A + beta*C.
  89 *>
  90 *>           Unchanged on exit.
  91 *> \endverbatim
  92 *>
  93 *> \param[in] N
  94 *> \verbatim
  95 *>          N is INTEGER
  96 *>           On entry, N specifies the order of the matrix C. N must be
  97 *>           at least zero.
  98 *>           Unchanged on exit.
  99 *> \endverbatim
 100 *>
 101 *> \param[in] K
 102 *> \verbatim
 103 *>          K is INTEGER
 104 *>           On entry with TRANS = 'N' or 'n', K specifies the number
 105 *>           of  columns of the matrix A, and on entry with TRANS = 'T'
 106 *>           or 't', K specifies the number of rows of the matrix A. K
 107 *>           must be at least zero.
 108 *>           Unchanged on exit.
 109 *> \endverbatim
 110 *>
 111 *> \param[in] ALPHA
 112 *> \verbatim
 113 *>          ALPHA is REAL
 114 *>           On entry, ALPHA specifies the scalar alpha.
 115 *>           Unchanged on exit.
 116 *> \endverbatim
 117 *>
 118 *> \param[in] A
 119 *> \verbatim
 120 *>          A is REAL array of DIMENSION (LDA,ka)
 121 *>           where KA
 122 *>           is K  when TRANS = 'N' or 'n', and is N otherwise. Before
 123 *>           entry with TRANS = 'N' or 'n', the leading N--by--K part of
 124 *>           the array A must contain the matrix A, otherwise the leading
 125 *>           K--by--N part of the array A must contain the matrix A.
 126 *>           Unchanged on exit.
 127 *> \endverbatim
 128 *>
 129 *> \param[in] LDA
 130 *> \verbatim
 131 *>          LDA is INTEGER
 132 *>           On entry, LDA specifies the first dimension of A as declared
 133 *>           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
 134 *>           then  LDA must be at least  max( 1, n ), otherwise  LDA must
 135 *>           be at least  max( 1, k ).
 136 *>           Unchanged on exit.
 137 *> \endverbatim
 138 *>
 139 *> \param[in] BETA
 140 *> \verbatim
 141 *>          BETA is REAL
 142 *>           On entry, BETA specifies the scalar beta.
 143 *>           Unchanged on exit.
 144 *> \endverbatim
 145 *>
 146 *> \param[in,out] C
 147 *> \verbatim
 148 *>          C is REAL array, dimension (NT)
 149 *>           NT = N*(N+1)/2. On entry, the symmetric matrix C in RFP
 150 *>           Format. RFP Format is described by TRANSR, UPLO and N.
 151 *> \endverbatim
 152 *
 153 *  Authors:
 154 *  ========
 155 *
 156 *> \author Univ. of Tennessee
 157 *> \author Univ. of California Berkeley
 158 *> \author Univ. of Colorado Denver
 159 *> \author NAG Ltd.
 160 *
 161 *> \date September 2012
 162 *
 163 *> \ingroup realOTHERcomputational
 164 *
 165 *  =====================================================================
 166       SUBROUTINE SSFRK( TRANSR, UPLO, TRANS, N, K, ALPHA, A, LDA, BETA,
 167      $                  C )
 168 *
 169 *  -- LAPACK computational routine (version 3.4.2) --
 170 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 171 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 172 *     September 2012
 173 *
 174 *     .. Scalar Arguments ..
 175       REAL               ALPHA, BETA
 176       INTEGER            K, LDA, N
 177       CHARACTER          TRANS, TRANSR, UPLO
 178 *     ..
 179 *     .. Array Arguments ..
 180       REAL               A( LDA, * ), C( * )
 181 *     ..
 182 *
 183 *  =====================================================================
 184 *
 185 *     .. Parameters ..
 186       REAL               ONE, ZERO
 187       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
 188 *     ..
 189 *     .. Local Scalars ..
 190       LOGICAL            LOWER, NORMALTRANSR, NISODD, NOTRANS
 191       INTEGER            INFO, NROWA, J, NK, N1, N2
 192 *     ..
 193 *     .. External Functions ..
 194       LOGICAL            LSAME
 195       EXTERNAL           LSAME
 196 *     ..
 197 *     .. External Subroutines ..
 198       EXTERNAL           SGEMM, SSYRK, XERBLA
 199 *     ..
 200 *     .. Intrinsic Functions ..
 201       INTRINSIC          MAX
 202 *     ..
 203 *     .. Executable Statements ..
 204 *
 205 *     Test the input parameters.
 206 *
 207       INFO = 0
 208       NORMALTRANSR = LSAME( TRANSR, 'N' )
 209       LOWER = LSAME( UPLO, 'L' )
 210       NOTRANS = LSAME( TRANS, 'N' )
 211 *
 212       IF( NOTRANS ) THEN
 213          NROWA = N
 214       ELSE
 215          NROWA = K
 216       END IF
 217 *
 218       IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN
 219          INFO = -1
 220       ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
 221          INFO = -2
 222       ELSE IF( .NOT.NOTRANS .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
 223          INFO = -3
 224       ELSE IF( N.LT.0 ) THEN
 225          INFO = -4
 226       ELSE IF( K.LT.0 ) THEN
 227          INFO = -5
 228       ELSE IF( LDA.LT.MAX( 1, NROWA ) ) THEN
 229          INFO = -8
 230       END IF
 231       IF( INFO.NE.0 ) THEN
 232          CALL XERBLA( 'SSFRK ', -INFO )
 233          RETURN
 234       END IF
 235 *
 236 *     Quick return if possible.
 237 *
 238 *     The quick return case: ((ALPHA.EQ.0).AND.(BETA.NE.ZERO)) is not
 239 *     done (it is in SSYRK for example) and left in the general case.
 240 *
 241       IF( ( N.EQ.0 ) .OR. ( ( ( ALPHA.EQ.ZERO ) .OR. ( K.EQ.0 ) ) .AND.
 242      $    ( BETA.EQ.ONE ) ) )RETURN
 243 *
 244       IF( ( ALPHA.EQ.ZERO ) .AND. ( BETA.EQ.ZERO ) ) THEN
 245          DO J = 1, ( ( N*( N+1 ) ) / 2 )
 246             C( J ) = ZERO
 247          END DO
 248          RETURN
 249       END IF
 250 *
 251 *     C is N-by-N.
 252 *     If N is odd, set NISODD = .TRUE., and N1 and N2.
 253 *     If N is even, NISODD = .FALSE., and NK.
 254 *
 255       IF( MOD( N, 2 ).EQ.0 ) THEN
 256          NISODD = .FALSE.
 257          NK = N / 2
 258       ELSE
 259          NISODD = .TRUE.
 260          IF( LOWER ) THEN
 261             N2 = N / 2
 262             N1 = N - N2
 263          ELSE
 264             N1 = N / 2
 265             N2 = N - N1
 266          END IF
 267       END IF
 268 *
 269       IF( NISODD ) THEN
 270 *
 271 *        N is odd
 272 *
 273          IF( NORMALTRANSR ) THEN
 274 *
 275 *           N is odd and TRANSR = 'N'
 276 *
 277             IF( LOWER ) THEN
 278 *
 279 *              N is odd, TRANSR = 'N', and UPLO = 'L'
 280 *
 281                IF( NOTRANS ) THEN
 282 *
 283 *                 N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
 284 *
 285                   CALL SSYRK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
 286      $                        BETA, C( 1 ), N )
 287                   CALL SSYRK( 'U', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
 288      $                        BETA, C( N+1 ), N )
 289                   CALL SGEMM( 'N', 'T', N2, N1, K, ALPHA, A( N1+1, 1 ),
 290      $                        LDA, A( 1, 1 ), LDA, BETA, C( N1+1 ), N )
 291 *
 292                ELSE
 293 *
 294 *                 N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'T'
 295 *
 296                   CALL SSYRK( 'L', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
 297      $                        BETA, C( 1 ), N )
 298                   CALL SSYRK( 'U', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
 299      $                        BETA, C( N+1 ), N )
 300                   CALL SGEMM( 'T', 'N', N2, N1, K, ALPHA, A( 1, N1+1 ),
 301      $                        LDA, A( 1, 1 ), LDA, BETA, C( N1+1 ), N )
 302 *
 303                END IF
 304 *
 305             ELSE
 306 *
 307 *              N is odd, TRANSR = 'N', and UPLO = 'U'
 308 *
 309                IF( NOTRANS ) THEN
 310 *
 311 *                 N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
 312 *
 313                   CALL SSYRK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
 314      $                        BETA, C( N2+1 ), N )
 315                   CALL SSYRK( 'U', 'N', N2, K, ALPHA, A( N2, 1 ), LDA,
 316      $                        BETA, C( N1+1 ), N )
 317                   CALL SGEMM( 'N', 'T', N1, N2, K, ALPHA, A( 1, 1 ),
 318      $                        LDA, A( N2, 1 ), LDA, BETA, C( 1 ), N )
 319 *
 320                ELSE
 321 *
 322 *                 N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'T'
 323 *
 324                   CALL SSYRK( 'L', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
 325      $                        BETA, C( N2+1 ), N )
 326                   CALL SSYRK( 'U', 'T', N2, K, ALPHA, A( 1, N2 ), LDA,
 327      $                        BETA, C( N1+1 ), N )
 328                   CALL SGEMM( 'T', 'N', N1, N2, K, ALPHA, A( 1, 1 ),
 329      $                        LDA, A( 1, N2 ), LDA, BETA, C( 1 ), N )
 330 *
 331                END IF
 332 *
 333             END IF
 334 *
 335          ELSE
 336 *
 337 *           N is odd, and TRANSR = 'T'
 338 *
 339             IF( LOWER ) THEN
 340 *
 341 *              N is odd, TRANSR = 'T', and UPLO = 'L'
 342 *
 343                IF( NOTRANS ) THEN
 344 *
 345 *                 N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
 346 *
 347                   CALL SSYRK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
 348      $                        BETA, C( 1 ), N1 )
 349                   CALL SSYRK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
 350      $                        BETA, C( 2 ), N1 )
 351                   CALL SGEMM( 'N', 'T', N1, N2, K, ALPHA, A( 1, 1 ),
 352      $                        LDA, A( N1+1, 1 ), LDA, BETA,
 353      $                        C( N1*N1+1 ), N1 )
 354 *
 355                ELSE
 356 *
 357 *                 N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
 358 *
 359                   CALL SSYRK( 'U', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
 360      $                        BETA, C( 1 ), N1 )
 361                   CALL SSYRK( 'L', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
 362      $                        BETA, C( 2 ), N1 )
 363                   CALL SGEMM( 'T', 'N', N1, N2, K, ALPHA, A( 1, 1 ),
 364      $                        LDA, A( 1, N1+1 ), LDA, BETA,
 365      $                        C( N1*N1+1 ), N1 )
 366 *
 367                END IF
 368 *
 369             ELSE
 370 *
 371 *              N is odd, TRANSR = 'T', and UPLO = 'U'
 372 *
 373                IF( NOTRANS ) THEN
 374 *
 375 *                 N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
 376 *
 377                   CALL SSYRK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
 378      $                        BETA, C( N2*N2+1 ), N2 )
 379                   CALL SSYRK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
 380      $                        BETA, C( N1*N2+1 ), N2 )
 381                   CALL SGEMM( 'N', 'T', N2, N1, K, ALPHA, A( N1+1, 1 ),
 382      $                        LDA, A( 1, 1 ), LDA, BETA, C( 1 ), N2 )
 383 *
 384                ELSE
 385 *
 386 *                 N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
 387 *
 388                   CALL SSYRK( 'U', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
 389      $                        BETA, C( N2*N2+1 ), N2 )
 390                   CALL SSYRK( 'L', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
 391      $                        BETA, C( N1*N2+1 ), N2 )
 392                   CALL SGEMM( 'T', 'N', N2, N1, K, ALPHA, A( 1, N1+1 ),
 393      $                        LDA, A( 1, 1 ), LDA, BETA, C( 1 ), N2 )
 394 *
 395                END IF
 396 *
 397             END IF
 398 *
 399          END IF
 400 *
 401       ELSE
 402 *
 403 *        N is even
 404 *
 405          IF( NORMALTRANSR ) THEN
 406 *
 407 *           N is even and TRANSR = 'N'
 408 *
 409             IF( LOWER ) THEN
 410 *
 411 *              N is even, TRANSR = 'N', and UPLO = 'L'
 412 *
 413                IF( NOTRANS ) THEN
 414 *
 415 *                 N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
 416 *
 417                   CALL SSYRK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
 418      $                        BETA, C( 2 ), N+1 )
 419                   CALL SSYRK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
 420      $                        BETA, C( 1 ), N+1 )
 421                   CALL SGEMM( 'N', 'T', NK, NK, K, ALPHA, A( NK+1, 1 ),
 422      $                        LDA, A( 1, 1 ), LDA, BETA, C( NK+2 ),
 423      $                        N+1 )
 424 *
 425                ELSE
 426 *
 427 *                 N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'T'
 428 *
 429                   CALL SSYRK( 'L', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
 430      $                        BETA, C( 2 ), N+1 )
 431                   CALL SSYRK( 'U', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
 432      $                        BETA, C( 1 ), N+1 )
 433                   CALL SGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, NK+1 ),
 434      $                        LDA, A( 1, 1 ), LDA, BETA, C( NK+2 ),
 435      $                        N+1 )
 436 *
 437                END IF
 438 *
 439             ELSE
 440 *
 441 *              N is even, TRANSR = 'N', and UPLO = 'U'
 442 *
 443                IF( NOTRANS ) THEN
 444 *
 445 *                 N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
 446 *
 447                   CALL SSYRK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
 448      $                        BETA, C( NK+2 ), N+1 )
 449                   CALL SSYRK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
 450      $                        BETA, C( NK+1 ), N+1 )
 451                   CALL SGEMM( 'N', 'T', NK, NK, K, ALPHA, A( 1, 1 ),
 452      $                        LDA, A( NK+1, 1 ), LDA, BETA, C( 1 ),
 453      $                        N+1 )
 454 *
 455                ELSE
 456 *
 457 *                 N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'T'
 458 *
 459                   CALL SSYRK( 'L', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
 460      $                        BETA, C( NK+2 ), N+1 )
 461                   CALL SSYRK( 'U', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
 462      $                        BETA, C( NK+1 ), N+1 )
 463                   CALL SGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, 1 ),
 464      $                        LDA, A( 1, NK+1 ), LDA, BETA, C( 1 ),
 465      $                        N+1 )
 466 *
 467                END IF
 468 *
 469             END IF
 470 *
 471          ELSE
 472 *
 473 *           N is even, and TRANSR = 'T'
 474 *
 475             IF( LOWER ) THEN
 476 *
 477 *              N is even, TRANSR = 'T', and UPLO = 'L'
 478 *
 479                IF( NOTRANS ) THEN
 480 *
 481 *                 N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
 482 *
 483                   CALL SSYRK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
 484      $                        BETA, C( NK+1 ), NK )
 485                   CALL SSYRK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
 486      $                        BETA, C( 1 ), NK )
 487                   CALL SGEMM( 'N', 'T', NK, NK, K, ALPHA, A( 1, 1 ),
 488      $                        LDA, A( NK+1, 1 ), LDA, BETA,
 489      $                        C( ( ( NK+1 )*NK )+1 ), NK )
 490 *
 491                ELSE
 492 *
 493 *                 N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
 494 *
 495                   CALL SSYRK( 'U', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
 496      $                        BETA, C( NK+1 ), NK )
 497                   CALL SSYRK( 'L', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
 498      $                        BETA, C( 1 ), NK )
 499                   CALL SGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, 1 ),
 500      $                        LDA, A( 1, NK+1 ), LDA, BETA,
 501      $                        C( ( ( NK+1 )*NK )+1 ), NK )
 502 *
 503                END IF
 504 *
 505             ELSE
 506 *
 507 *              N is even, TRANSR = 'T', and UPLO = 'U'
 508 *
 509                IF( NOTRANS ) THEN
 510 *
 511 *                 N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
 512 *
 513                   CALL SSYRK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
 514      $                        BETA, C( NK*( NK+1 )+1 ), NK )
 515                   CALL SSYRK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
 516      $                        BETA, C( NK*NK+1 ), NK )
 517                   CALL SGEMM( 'N', 'T', NK, NK, K, ALPHA, A( NK+1, 1 ),
 518      $                        LDA, A( 1, 1 ), LDA, BETA, C( 1 ), NK )
 519 *
 520                ELSE
 521 *
 522 *                 N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
 523 *
 524                   CALL SSYRK( 'U', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
 525      $                        BETA, C( NK*( NK+1 )+1 ), NK )
 526                   CALL SSYRK( 'L', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
 527      $                        BETA, C( NK*NK+1 ), NK )
 528                   CALL SGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, NK+1 ),
 529      $                        LDA, A( 1, 1 ), LDA, BETA, C( 1 ), NK )
 530 *
 531                END IF
 532 *
 533             END IF
 534 *
 535          END IF
 536 *
 537       END IF
 538 *
 539       RETURN
 540 *
 541 *     End of SSFRK
 542 *
 543       END