Mdcxmatinv.cxx

Go to the documentation of this file.

#include <math.h>
#include <iostream>
#include "MdcxReco/Mdcxmatinv.h"
using std::cout;
using std::endl;
 
extern int Mdcxmatinv(double *array, int *norder, double *det){
    /* System generated locals */
    int i__3;
    double d__1;
 
    /* Local variables */
    const int nmax = 10;
//    cout << "norder in Mdcxmatinv = " << *norder << endl;
    if (*norder > nmax){
     cout << "In Mdcxmatinv, norder ( = " << *norder << " ) > nmax ( = "
          << nmax << " ); error" << endl; return 1000;
    }
    static double amax, save;
    static int i, j, k, l, ik[nmax], jk[nmax];
 
    /* Parameter adjustments */
    array -= (nmax+1);
 
    /* Function Body */
    *det = (double)1.;
    for (k = 1; k <= *norder; ++k) {
 
/*       FIND LARGEST ELEMENT ARRAY(I, J) IN REST OF MATRIX */
 
    amax = (double)0.;
L21:
    for (i = k; i <= *norder; ++i) {
        for (j = k; j <= *norder; ++j) {
        d__1 = array[i + j * nmax]; 
        if ((fabs(amax)-fabs(d__1)) <= 0.) {
            amax = array[i + j * nmax];
            ik[k - 1] = i;
            jk[k - 1] = j;
        }
        }
    }
 
/*       INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN ARRAY(K, K) */
 
    if (amax == 0.) {*det = (double)0.; return 1001;}
 
    i = ik[k - 1];
    if ((i__3 = i - k) < 0) {
        goto L21;
    } else if (i__3 == 0) {
        goto L51;
    } else {
        goto L43;
    }
L43:
    for (j = 1; j <= *norder; ++j) {
        save = array[k + j * nmax];
        array[k + j * nmax] = array[i + j * nmax];
        array[i + j * nmax] = -save;
    }
L51:
    j = jk[k - 1];
    if ((i__3 = j - k) < 0) {
        goto L21;
    } else if (i__3 == 0) {
        goto L61;
    } else {
        goto L53;
    }
L53:
    for (i = 1; i <= *norder; ++i) {
        save = array[i + k * nmax];
        array[i + k * nmax] = array[i + j * nmax];
        array[i + j * nmax] = -save;
    }
 
/*       ACCUMULATE ELEMENTS OF INVERSE MATRIX */
 
L61:
    for (i = 1; i <= *norder; ++i) {
        if (i - k != 0) {
            array[i + k * nmax] = -array[i + k * nmax] / amax;
        }   
    }
    for (i = 1; i <= *norder; ++i) {
        for (j = 1; j <= *norder; ++j) {
        if (i - k != 0) {
            goto L74;
        } else {
            goto L80;
        }
L74:
        if (j - k != 0) {
            goto L75;
        } else {
            goto L80;
        }
L75:
        array[i+j*nmax] += array[i+k*nmax] * array[k+j*nmax];
L80:
        ;
        }
    }
    for (j = 1; j <= *norder; ++j) {
        if (j - k != 0) {
            array[k + j * nmax] /= amax;
        }
    }
    array[k + k * nmax] = (double)1. / amax;
    *det *= amax;
    }
 
/*       RESTORE ORDERING OF MATRIX */
 
    for (l = 1; l <= *norder; ++l) {
    k = *norder - l + 1;
    j = ik[k - 1];
    if (j - k <= 0) {
        goto L111;
    } else {
        goto L105;
    }
L105:
    for (i = 1; i <= *norder; ++i) {
        save = array[i + k * nmax];
        array[i + k * nmax] = -array[i + j * nmax];
        array[i + j * nmax] = save;
    }
L111:
    i = jk[k - 1];
    if (i - k <= 0) {
        goto L130;
    } else {
        goto L113;
    }
L113:
    for (j = 1; j <= *norder; ++j) {
        save = array[k + j * nmax];
        array[k + j * nmax] = -array[i + j * nmax];
        array[i + j * nmax] = save;
    }
L130:
    ;
    }
    return 0;
} /* Mdcxmatinv */