00001 #include <math.h>
00002 #include <iostream>
00003 #include "MdcxReco/Mdcxmatinv.h"
00004 using std::cout;
00005 using std::endl;
00006
00007 extern int Mdcxmatinv(double *array, int *norder, double *det){
00008
00009 int i__3;
00010 double d__1;
00011
00012
00013 const int nmax = 10;
00014
00015 if (*norder > nmax){
00016 cout << "In Mdcxmatinv, norder ( = " << *norder << " ) > nmax ( = "
00017 << nmax << " ); error" << endl; return 1000;
00018 }
00019 static double amax, save;
00020 static int i, j, k, l, ik[nmax], jk[nmax];
00021
00022
00023 array -= (nmax+1);
00024
00025
00026 *det = (double)1.;
00027 for (k = 1; k <= *norder; ++k) {
00028
00029
00030
00031 amax = (double)0.;
00032 L21:
00033 for (i = k; i <= *norder; ++i) {
00034 for (j = k; j <= *norder; ++j) {
00035 d__1 = array[i + j * nmax];
00036 if ((fabs(amax)-fabs(d__1)) <= 0.) {
00037 amax = array[i + j * nmax];
00038 ik[k - 1] = i;
00039 jk[k - 1] = j;
00040 }
00041 }
00042 }
00043
00044
00045
00046 if (amax == 0.) {*det = (double)0.; return 1001;}
00047
00048 i = ik[k - 1];
00049 if ((i__3 = i - k) < 0) {
00050 goto L21;
00051 } else if (i__3 == 0) {
00052 goto L51;
00053 } else {
00054 goto L43;
00055 }
00056 L43:
00057 for (j = 1; j <= *norder; ++j) {
00058 save = array[k + j * nmax];
00059 array[k + j * nmax] = array[i + j * nmax];
00060 array[i + j * nmax] = -save;
00061 }
00062 L51:
00063 j = jk[k - 1];
00064 if ((i__3 = j - k) < 0) {
00065 goto L21;
00066 } else if (i__3 == 0) {
00067 goto L61;
00068 } else {
00069 goto L53;
00070 }
00071 L53:
00072 for (i = 1; i <= *norder; ++i) {
00073 save = array[i + k * nmax];
00074 array[i + k * nmax] = array[i + j * nmax];
00075 array[i + j * nmax] = -save;
00076 }
00077
00078
00079
00080 L61:
00081 for (i = 1; i <= *norder; ++i) {
00082 if (i - k != 0) {
00083 array[i + k * nmax] = -array[i + k * nmax] / amax;
00084 }
00085 }
00086 for (i = 1; i <= *norder; ++i) {
00087 for (j = 1; j <= *norder; ++j) {
00088 if (i - k != 0) {
00089 goto L74;
00090 } else {
00091 goto L80;
00092 }
00093 L74:
00094 if (j - k != 0) {
00095 goto L75;
00096 } else {
00097 goto L80;
00098 }
00099 L75:
00100 array[i+j*nmax] += array[i+k*nmax] * array[k+j*nmax];
00101 L80:
00102 ;
00103 }
00104 }
00105 for (j = 1; j <= *norder; ++j) {
00106 if (j - k != 0) {
00107 array[k + j * nmax] /= amax;
00108 }
00109 }
00110 array[k + k * nmax] = (double)1. / amax;
00111 *det *= amax;
00112 }
00113
00114
00115
00116 for (l = 1; l <= *norder; ++l) {
00117 k = *norder - l + 1;
00118 j = ik[k - 1];
00119 if (j - k <= 0) {
00120 goto L111;
00121 } else {
00122 goto L105;
00123 }
00124 L105:
00125 for (i = 1; i <= *norder; ++i) {
00126 save = array[i + k * nmax];
00127 array[i + k * nmax] = -array[i + j * nmax];
00128 array[i + j * nmax] = save;
00129 }
00130 L111:
00131 i = jk[k - 1];
00132 if (i - k <= 0) {
00133 goto L130;
00134 } else {
00135 goto L113;
00136 }
00137 L113:
00138 for (j = 1; j <= *norder; ++j) {
00139 save = array[k + j * nmax];
00140 array[k + j * nmax] = -array[i + j * nmax];
00141 array[i + j * nmax] = save;
00142 }
00143 L130:
00144 ;
00145 }
00146 return 0;
00147 }