1

In numerical recipes, on page 39 (page 4 of the pdf) the following algorithm has been suggested for finding a pivot:

void gaussj(float **a, int n, float **b, int m)
/* Linear equation solution by Gauss-Jordan elimination,
   equation (2.1.1) above. a[1..n][1..n] is the input matrix.
   b[1..n][1..m] is input containing the m right-hand side vectors.
   On output, a is replaced by its matrix inverse, and b is
   replaced by the corresponding set of solution vectors.  */
{
   int *indxc,*indxr,*ipiv;
   int i,icol,irow,j,k,l,ll;
   float big,dum,pivinv,temp;
   indxc=ivector(1,n); indxr=ivector(1,n); ipiv=ivector(1,n);
   /* The integer arrays ipiv, indxr, and indxc are used for
      bookkeeping on the pivoting. */
   for (j=1;j<=n;j++) ipiv[j]=0;
   // This is the main loop over the columns to be reduced.
   for  (i=1;i<=n;i++) {
       big=0.0;
       // The following is the outer loop of the
       // search for a pivot element.
       for (j=1;j<=n;j++)
           if (ipiv[j] != 1)
               for (k=1;k<=n;k++) {
                   if (ipiv[k] == 0) {
                       if (abs(a[j][k]) >= big) {
                           big=abs(a[j][k]);
                               irow=j;
                               icol=k;

                       }
                   }
               }
       ++(ipiv[icol]);

       // ... THE REST OF THE CODE ... (complete code in the link)

I don't understand why at some point we have ipiv[j] != 1 and two lines later ipiv[k] == 0. And then at end an element of the vector is incremented. What is the semantic of this vector?

Just to expand on this, is there any reason that one line compares an element of ipiv to 1 and two lines later it is compared against 0. Also the last line increments an element instead of simply setting it to 1. So is this just a lack of consistency and style in writing the code, or an element of ipiv can take values other than 0 and 1?

Edit: I don't know how to activate syntax-highlighting for the code. If you know how, please let me know.

igor
  • 473
  • I don't think I understand the layout of the braces in this code: where is the open brace for the j loop? I also think there are several missing close braces... – Ian Feb 21 '16 at 02:34
  • @Ian The complete code is in the link. I fixed an indentation issue. The "missing" braces, I think, are related to C syntax where if the body of a loop or if has only one statement, you can drop the braces. – igor Feb 21 '16 at 02:49
  • I guess you loop over columns i, then loop over rows j, if you already have a pivot for the row j then you skip, otherwise you loop over columns to find a pivot, jut down what you found, and then update that you've found a pivot for the column in question. Or something like that, this is not really written for readability as far as I can tell. – Ian Feb 21 '16 at 03:05
  • Yes, it is clearer in the full code: irow,icol is the position of the pivot being found in a given iteration of the loop, and ipiv tells you whether a pivot has already been found for a given row j. – Ian Feb 21 '16 at 03:07
  • @Ian do you think that an element of ipiv can either be zero or one? Is there any reason that one line compares an element of ipiv to 1 and two lines later it is compared against 0. Also the last line increments the value instead of setting it to 1. So is this just a lack of consistency and style in writing the code, or an element of ipiv can take values other than 0 and 1? – igor Feb 21 '16 at 03:32
  • Assuming it's only used in an if, incrementing it again does no real harm, it's still a "true". Anyway, two lines later you are looking at a generally different element of ipiv, indexed by k inside the new k loop. – Ian Feb 21 '16 at 03:59

0 Answers0