public static void choldc_f77(int n, double a[][], double diagmx,
                                double tol, double addmax[]) {

/*

Here is a copy of the choldc FORTRAN documentation:

      SUBROUTINE CHOLDC(NR,N,A,DIAGMX,TOL,ADDMAX)

C
C PURPOSE
C -------
C FIND THE PERTURBED L(L-TRANSPOSE) [WRITTEN LL+] DECOMPOSITION
C OF A+D, WHERE D IS A NON-NEGATIVE DIAGONAL MATRIX ADDED TO A IF
C NECESSARY TO ALLOW THE CHOLESKY DECOMPOSITION TO CONTINUE.
C
C PARAMETERS
C ----------
C NR           --> ROW DIMENSION OF MATRIX
C N            --> DIMENSION OF PROBLEM
C A(N,N)      <--> ON ENTRY: MATRIX FOR WHICH TO FIND PERTURBED
C                       CHOLESKY DECOMPOSITION
C                  ON EXIT:  CONTAINS L OF LL+ DECOMPOSITION
C                  IN LOWER TRIANGULAR PART AND DIAGONAL OF "A"
C DIAGMX       --> MAXIMUM DIAGONAL ELEMENT OF "A"
C TOL          --> TOLERANCE
C ADDMAX      <--  MAXIMUM AMOUNT IMPLICITLY ADDED TO DIAGONAL OF "A"
C                  IN FORMING THE CHOLESKY DECOMPOSITION OF A+D
C INTERNAL VARIABLES
C ------------------
C AMINL    SMALLEST ELEMENT ALLOWED ON DIAGONAL OF L
C AMNLSQ   =AMINL**2
C OFFMAX   MAXIMUM OFF-DIAGONAL ELEMENT IN COLUMN OF A
C
C
C DESCRIPTION
C -----------
C THE NORMAL CHOLESKY DECOMPOSITION IS PERFORMED.  HOWEVER, IF AT ANY
C POINT THE ALGORITHM WOULD ATTEMPT TO SET L(I,I)=SQRT(TEMP)
C WITH TEMP < TOL*DIAGMX, THEN L(I,I) IS SET TO SQRT(TOL*DIAGMX)
C INSTEAD.  THIS IS EQUIVALENT TO ADDING TOL*DIAGMX-TEMP TO A(I,I)
C
C

*/

      int i,j,jm1,jp1,k;
      double aminl,amnlsq,offmax,sum,temp;

      addmax[1] = 0.0;
      aminl = Math.sqrt(diagmx*tol);
      amnlsq = aminl*aminl;

// FORM COLUMN J OF L

      for (j = 1; j <= n; j++) {

// FIND DIAGONAL ELEMENTS OF L

         sum = 0.0;
         jm1 = j - 1;
         jp1 = j + 1;

         for (k = 1; k <= jm1; k++) {

            sum += a[j][k]*a[j][k];

         }

         temp = a[j][j] - sum;

         if (temp >= amnlsq) {

            a[j][j] = Math.sqrt(temp);

         } else {

// FIND MAXIMUM OFF-DIAGONAL ELEMENT IN COLUMN

            offmax = 0.0;

            for (i = jp1; i <= n; i++) {

               if (Math.abs(a[i][j]) > offmax) offmax = Math.abs(a[i][j]);

            }

            if (offmax <= amnlsq) offmax = amnlsq;

// ADD TO DIAGONAL ELEMENT  TO ALLOW CHOLESKY DECOMPOSITION TO CONTINUE

            a[j][j] = Math.sqrt(offmax);
            addmax[1] = Math.max(addmax[1],offmax-temp);

         }

// FIND I,J ELEMENT OF LOWER TRIANGULAR MATRIX

         for (i = jp1; i <= n; i++) {

            sum = 0.0;

            for (k = 1; k <= jm1; k++) {

               sum += a[i][k]*a[j][k];

            }

            a[i][j] = (a[i][j] - sum)/a[j][j];

         }

      }

      return;

   }