Here is an example from the IMSL library of mathematical functions. The function DGEFS (which solves a set of linear equations) has 105 lines of comments and 60 lines of program instructions.

SUBROUTINE DGEFS (A, LDA, N, V, ITASK, IND, WORK, IWORK) C***BEGIN PROLOGUE DGEFS C***PURPOSE Solve a general system of linear equations. C***LIBRARY SLATEC C***CATEGORY D2A1 C***TYPE DOUBLE PRECISION (SGEFS-S, DGEFS-D, CGEFS-C) C***KEYWORDS COMPLEX LINEAR EQUATIONS, GENERAL MATRIX, C GENERAL SYSTEM OF LINEAR EQUATIONS C***AUTHOR Voorhees, E. A., (LANL) C***DESCRIPTION C C Subroutine DGEFS solves a general NxN system of double C precision linear equations using LINPACK subroutines DGECO C and DGESL. That is, if A is an NxN double precision matrix C and if X and B are double precision N-vectors, then DGEFS C solves the equation C C A*X=B. C C The matrix A is first factored into upper and lower tri- C angular matrices U and L using partial pivoting. These C factors and the pivoting information are used to find the C solution vector X. An approximate condition number is C calculated to provide a rough estimate of the number of C digits of accuracy in the computed solution. C C If the equation A*X=B is to be solved for more than one vector C B, the factoring of A does not need to be performed again and C the option to only solve (ITASK.GT.1) will be faster for C the succeeding solutions. In this case, the contents of A, C LDA, N and IWORK must not have been altered by the user follow- C ing factorization (ITASK=1). IND will not be changed by DGEFS C in this case. C C Argument Description *** C C A DOUBLE PRECISION(LDA,N) C on entry, the doubly subscripted array with dimension C (LDA,N) which contains the coefficient matrix. C on return, an upper triangular matrix U and the C multipliers necessary to construct a matrix L C so that A=L*U. C LDA INTEGER C the leading dimension of the array A. LDA must be great- C er than or equal to N. (terminal error message IND=-1) C N INTEGER C the order of the matrix A. The first N elements of C the array A are the elements of the first column of C the matrix A. N must be greater than or equal to 1. C (terminal error message IND=-2) C V DOUBLE PRECISION(N) C on entry, the singly subscripted array(vector) of di- C mension N which contains the right hand side B of a C system of simultaneous linear equations A*X=B. C on return, V contains the solution vector, X . C ITASK INTEGER C If ITASK=1, the matrix A is factored and then the C linear equation is solved. C If ITASK .GT. 1, the equation is solved using the existing C factored matrix A and IWORK. C If ITASK .LT. 1, then terminal error message IND=-3 is C printed. C IND INTEGER C GT. 0 IND is a rough estimate of the number of digits C of accuracy in the solution, X. C LT. 0 see error message corresponding to IND below. C WORK DOUBLE PRECISION(N) C a singly subscripted array of dimension at least N. C IWORK INTEGER(N) C a singly subscripted array of dimension at least N. C C Error Messages Printed *** C C IND=-1 terminal N is greater than LDA. C IND=-2 terminal N is less than 1. C IND=-3 terminal ITASK is less than 1. C IND=-4 terminal The matrix A is computationally singular. C A solution has not been computed. C IND=-10 warning The solution has no apparent significance. C The solution may be inaccurate or the matrix C A may be poorly scaled. C C Note- The above terminal(*fatal*) error messages are C designed to be handled by XERMSG in which C LEVEL=1 (recoverable) and IFLAG=2 . LEVEL=0 C for warning error messages from XERMSG. Unless C the user provides otherwise, an error message C will be printed followed by an abort. C C***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W. C Stewart, LINPACK Users' Guide, SIAM, 1979. C***ROUTINES CALLED D1MACH, DGECO, DGESL, XERMSG C***REVISION HISTORY (YYMMDD) C 800326 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890831 Modified array declarations. (WRB) C 890831 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) C 900510 Convert XERRWV calls to XERMSG calls. (RWC) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE DGEFS C the Fortran code starts here ....

An example from a C program I wrote: about 53 lines of comments and description, followed by about 60 lines of program instructions.

/************************************************************************ ** * ROUTINE: atMatchLists * * DESCRIPTION: * Given 2 lists of s_star structures, * which have ALREADY been transformed so that the "x" * and "y" coordinates of each list are close to each other * (i.e. matching items from each list have very similar "x" and "y") * this routine attempts to find all instances of matching items * from the 2 lists. * * We consider a "match" to be the closest coincidence of centers * which are within "radius" pixels of each other. * * Use a slow, but sure, algorithm. * * We will match objects from A --> B. It is possible to have several * As that match to the same B: * * A1 -> B5 and A2 -> B5 * * This function finds such multiple-match items and deletes all but * the closest of the matches. * * place the elems of A that are matches into output list J * B that are matches into output list K * A that are not matches into output list L * B that are not matches into output list M * * * RETURN: * SH_SUCCESS if all goes well * SH_GENERIC_ERROR if an error occurs * * */ int atMatchLists ( int numA, /* I: number of stars in list A */ s_star *listA, /* I: first input list of items to be matched */ int numB, /* I: number of stars in list B */ s_star *listB, /* I: second list of items to be matched */ double radius, /* I: maximum radius for items to be a match */ char *basename /* I: base of filenames used to store the */ /* output; extension indicates which */ /* .mtA items from A that matched */ /* .mtB items from B that matched */ /* .unA items from A that didn't match */ /* .unB items from A that didn't match */ ) code starts here ...

I worked on a project which involved writing thousands of lines of code, broken up into over 50 files. In one of those files was this function, which is supposed to invert a matrix. It took quite a while for me to figure out what was going on, and how to use the routine. Grrrrr.

/************************************************************************** * * * inverts a symmetric routine * * * **************************************************************************/ /* Return 1 for success, 0 for failure */ #include/* Prototype definition */ #include "mtLstsq.h" int syminv(double *a, double *b, int dimension, double *pvrow, double *pvcol, int *use) code starts here ....