• Choose a category
• All
• Classic-Fiction

By Root 5207 0

that a computer has limitations. As I said earlier in this chapter, computers have limited speed, sometimes operate better on floating-point numbers or integer numbers, and have finite amounts of memory. Beautiful code must consider these limitations of reality. Quite often, people writing code assume that memory is infinite, computer speed is infinite, and so on. This is not beautiful code; it's arrogant code. Beautiful code is frugal about things like memory use and reuses memory whenever possible.

For example, let's look at the LU decomposition subroutine SGBTRF, which is in the second level of subroutines. To save space, I removed the initial comments in the header and other excepts that I do not directly discuss (you can view the complete subroutine at http://www.netlib.org/lapack/explore-html/sgbtrf.f.html):

Code View: Scroll / Show All

SUBROUTINE SGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO )

*

* -- LAPACK routine (version 2.0) --

* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,

* Courant Institute, Argonne National Lab, and Rice University

* February 29, 1992

.

.

.

Initial comments, description of parameters

.

,

,

* Test the input parameters.

*

INFO = 0

IF( M.LT.0 ) THEN

.

.

.

Checking parameters

.

.

.

CALL XERBLA( 'SGBTRF', -INFO )

RETURN

END IF

*

* Quick return if possible

*

IF( M.EQ.0 .OR. N.EQ.0 )

$ RETURN

*

* Determine the block size for this environment

*

NB = ILAENV( 1, 'SGBTRF', ' ', M, N, KL, KU )

*

* The block size must not exceed the limit set by the size of the

* local arrays WORK13 and WORK31.

*

NB = MIN( NB, NBMAX )

*

IF( NB.LE.1 .OR. NB.GT.KL ) THEN

*

* Use unblocked code

*

CALL SGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO )

ELSE

*

* Use blocked code

*

* Zero the superdiagonal elements of the work array WORK13

*

DO 20 J = 1, NB

DO 10 I = 1, J - 1

WORK13( I, J ) = ZERO

10 CONTINUE

20 CONTINUE

*

* Zero the subdiagonal elements of the work array WORK31

*

DO 40 J = 1, NB

DO 30 I = J + 1, NB

WORK31( I, J ) = ZERO

30 CONTINUE

40 CONTINUE

*

* Gaussian elimination with partial pivoting

*

* Set fill-in elements in columns KU+2 to KV to zero

*

DO 60 J = KU + 2, MIN( KV, N )

DO 50 I = KV - J + 2, KL

AB( I, J ) = ZERO

50 CONTINUE

60 CONTINUE

*

* JU is the index of the last column affected by the current

* stage of the factorization

*

JU = 1

*

DO 180 J = 1, MIN( M, N ), NB

JB = MIN( NB, MIN( M, N )-J+1 )

*

* The active part of the matrix is partitioned

*

* A11 A12 A13

* A21 A22 A23

* A31 A32 A33

*

* Here A11, A21 and A31 denote the current block of JB columns

* which is about to be factorized. The number of rows in the

* partitioning are JB, I2, I3 respectively, and the numbers

* of columns are JB, J2, J3. The superdiagonal elements of A13

* and the subdiagonal elements of A31 lie outside the band.

*

I2 = MIN( KL-JB, M-J-JB+1 )

I3 = MIN( JB, M-J-KL+1 )

*

* J2 and J3 are computed after JU has been updated.

*

* Factorize the current block of JB columns

*

DO 80 JJ = J, J + JB - 1

*

* Set fill-in elements in column JJ+KV to zero

*

IF( JJ+KV.LE.N ) THEN

DO 70 I = 1, KL

AB( I, JJ+KV ) = ZERO

70 CONTINUE

END IF

*

* Find pivot and test for singularity. KM is the number of

* subdiagonal elements in the current column.

*

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

*

* Apply interchange to columns J to J+JB-1

*

IF( JP+JJ-1.LT.J+KL ) THEN

*

CALL SSWAP( JB, AB( KV+1+JJ-J, J ), LDAB-1,

$ AB( KV+JP+JJ-J, J ), LDAB-1 )

ELSE

*

* The interchange affects columns J to JJ-1 of A31

* which are stored in the work array WORK31

*

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,