sgbtrf.f.html

famous linear algebra library (LAPACK) ports to windows

<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span>      IF( M.EQ.0 .OR. N.EQ.0 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Determine the block size for this environment
</span><span class="comment">*</span><span class="comment">
</span>      NB = <a name="ILAENV.146"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="SGBTRF.146"></a><a href="sgbtrf.f.html#SGBTRF.1">SGBTRF</a>'</span>, <span class="string">' '</span>, M, N, KL, KU )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     The block size must not exceed the limit set by the size of the
</span><span class="comment">*</span><span class="comment">     local arrays WORK13 and WORK31.
</span><span class="comment">*</span><span class="comment">
</span>      NB = MIN( NB, NBMAX )
<span class="comment">*</span><span class="comment">
</span>      IF( NB.LE.1 .OR. NB.GT.KL ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Use unblocked code
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SGBTF2.157"></a><a href="sgbtf2.f.html#SGBTF2.1">SGBTF2</a>( M, N, KL, KU, AB, LDAB, IPIV, INFO )
      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Use blocked code
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Zero the superdiagonal elements of the work array WORK13
</span><span class="comment">*</span><span class="comment">
</span>         DO 20 J = 1, NB
            DO 10 I = 1, J - 1
               WORK13( I, J ) = ZERO
   10       CONTINUE
   20    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Zero the subdiagonal elements of the work array WORK31
</span><span class="comment">*</span><span class="comment">
</span>         DO 40 J = 1, NB
            DO 30 I = J + 1, NB
               WORK31( I, J ) = ZERO
   30       CONTINUE
   40    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Gaussian elimination with partial pivoting
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Set fill-in elements in columns KU+2 to KV to zero
</span><span class="comment">*</span><span class="comment">
</span>         DO 60 J = KU + 2, MIN( KV, N )
            DO 50 I = KV - J + 2, KL
               AB( I, J ) = ZERO
   50       CONTINUE
   60    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        JU is the index of the last column affected by the current
</span><span class="comment">*</span><span class="comment">        stage of the factorization
</span><span class="comment">*</span><span class="comment">
</span>         JU = 1
<span class="comment">*</span><span class="comment">
</span>         DO 180 J = 1, MIN( M, N ), NB
            JB = MIN( NB, MIN( M, N )-J+1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           The active part of the matrix is partitioned
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              A11   A12   A13
</span><span class="comment">*</span><span class="comment">              A21   A22   A23
</span><span class="comment">*</span><span class="comment">              A31   A32   A33
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Here A11, A21 and A31 denote the current block of JB columns
</span><span class="comment">*</span><span class="comment">           which is about to be factorized. The number of rows in the
</span><span class="comment">*</span><span class="comment">           partitioning are JB, I2, I3 respectively, and the numbers
</span><span class="comment">*</span><span class="comment">           of columns are JB, J2, J3. The superdiagonal elements of A13
</span><span class="comment">*</span><span class="comment">           and the subdiagonal elements of A31 lie outside the band.
</span><span class="comment">*</span><span class="comment">
</span>            I2 = MIN( KL-JB, M-J-JB+1 )
            I3 = MIN( JB, M-J-KL+1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           J2 and J3 are computed after JU has been updated.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Factorize the current block of JB columns
</span><span class="comment">*</span><span class="comment">
</span>            DO 80 JJ = J, J + JB - 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Set fill-in elements in column JJ+KV to zero
</span><span class="comment">*</span><span class="comment">
</span>               IF( JJ+KV.LE.N ) THEN
                  DO 70 I = 1, KL
                     AB( I, JJ+KV ) = ZERO
   70             CONTINUE
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Find pivot and test for singularity. KM is the number of
</span><span class="comment">*</span><span class="comment">              subdiagonal elements in the current column.
</span><span class="comment">*</span><span class="comment">
</span>               KM = MIN( KL, M-JJ )
               JP = ISAMAX( KM+1, AB( KV+1, JJ ), 1 )
               IPIV( JJ ) = JP + JJ - J
               IF( AB( KV+JP, JJ ).NE.ZERO ) THEN
                  JU = MAX( JU, MIN( JJ+KU+JP-1, N ) )
                  IF( JP.NE.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Apply interchange to columns J to J+JB-1
</span><span class="comment">*</span><span class="comment">
</span>                     IF( JP+JJ-1.LT.J+KL ) THEN
<span class="comment">*</span><span class="comment">
</span>                        CALL SSWAP( JB, AB( KV+1+JJ-J, J ), LDAB-1,
     $                              AB( KV+JP+JJ-J, J ), LDAB-1 )
                     ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                       The interchange affects columns J to JJ-1 of A31
</span><span class="comment">*</span><span class="comment">                       which are stored in the work array WORK31
</span><span class="comment">*</span><span class="comment">
</span>                        CALL SSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1,
     $                              WORK31( JP+JJ-J-KL, 1 ), LDWORK )
                        CALL SSWAP( J+JB-JJ, AB( KV+1, JJ ), LDAB-1,
     $                              AB( KV+JP, JJ ), LDAB-1 )
                     END IF
                  END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Compute multipliers
</span><span class="comment">*</span><span class="comment">
</span>                  CALL SSCAL( KM, ONE / AB( KV+1, JJ ), AB( KV+2, JJ ),
     $                        1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Update trailing submatrix within the band and within
</span><span class="comment">*</span><span class="comment">                 the current block. JM is the index of the last column
</span><span class="comment">*</span><span class="comment">                 which needs to be updated.
</span><span class="comment">*</span><span class="comment">
</span>                  JM = MIN( JU, J+JB-1 )
                  IF( JM.GT.JJ )
     $               CALL SGER( KM, JM-JJ, -ONE, AB( KV+2, JJ ), 1,
     $                          AB( KV, JJ+1 ), LDAB-1,
     $                          AB( KV+1, JJ+1 ), LDAB-1 )
               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 If pivot is zero, set INFO to the index of the pivot
</span><span class="comment">*</span><span class="comment">                 unless a zero pivot has already been found.
</span><span class="comment">*</span><span class="comment">
</span>                  IF( INFO.EQ.0 )
     $               INFO = JJ
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Copy current column of A31 into the work array WORK31
</span><span class="comment">*</span><span class="comment">
</span>               NW = MIN( JJ-J+1, I3 )
               IF( NW.GT.0 )
     $            CALL SCOPY( NW, AB( KV+KL+1-JJ+J, JJ ), 1,
     $                        WORK31( 1, JJ-J+1 ), 1 )
   80       CONTINUE
            IF( J+JB.LE.N ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Apply the row interchanges to the other blocks.
</span><span class="comment">*</span><span class="comment">
</span>               J2 = MIN( JU-J+1, KV ) - JB
               J3 = MAX( 0, JU-J-KV+1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Use <a name="SLASWP.290"></a><a href="slaswp.f.html#SLASWP.1">SLASWP</a> to apply the row interchanges to A12, A22, and
</span><span class="comment">*</span><span class="comment">              A32.
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="SLASWP.293"></a><a href="slaswp.f.html#SLASWP.1">SLASWP</a>( J2, AB( KV+1-JB, J+JB ), LDAB-1, 1, JB,