[cig-commits] r17194 - mc/3D/CitcomS/trunk/lib

becker at geodynamics.org becker at geodynamics.org
Wed Sep 15 03:40:50 PDT 2010


Author: becker
Date: 2010-09-15 03:40:50 -0700 (Wed, 15 Sep 2010)
New Revision: 17194

Modified:
   mc/3D/CitcomS/trunk/lib/General_matrix_functions.c
   mc/3D/CitcomS/trunk/lib/prototypes.h
Log:
Missed one routine.



Modified: mc/3D/CitcomS/trunk/lib/General_matrix_functions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/General_matrix_functions.c	2010-09-15 10:36:58 UTC (rev 17193)
+++ mc/3D/CitcomS/trunk/lib/General_matrix_functions.c	2010-09-15 10:40:50 UTC (rev 17194)
@@ -782,6 +782,36 @@
 
 
 
+double cofactor(A,i,j,n)
+     double A[4][4];
+     int i,j,n;
+
+{ int k,l,p,q;
+  double determinant();
+
+  double B[4][4]; /* because of recursive behaviour of det/cofac, need to use
+			       new copy of B at each 'n' level of this routine */
+
+  if (n>3) printf("Error, no cofactors for matrix more than 3x3\n");
+
+  p=q=1;
+
+  for(k=1;k<=n;k++)    {
+     if(k==i) continue;
+     for(l=1;l<=n;l++)      {
+	   if (l==j) continue;
+           B[p][q]=A[k][l];
+	   q++ ;
+	   }
+     q=1;p++;
+     }
+
+
+  return(epsilon[i][j]*determinant(B,n-1));
+
+
+}
+
 long double lg_pow(long double a, int n)
 {
     /* compute the value of "a" raised to the power of "n" */

Modified: mc/3D/CitcomS/trunk/lib/prototypes.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/prototypes.h	2010-09-15 10:36:58 UTC (rev 17193)
+++ mc/3D/CitcomS/trunk/lib/prototypes.h	2010-09-15 10:40:50 UTC (rev 17194)
@@ -130,6 +130,7 @@
 double conj_grad(struct All_variables *, double **, double **, double, int *, int);
 void element_gauss_seidel(struct All_variables *, double **, double **, double **, double, int *, int, int);
 void gauss_seidel(struct All_variables *, double **, double **, double **, double, int *, int, int);
+double cofactor(double [4][4], int, int, int);
 
 double determinant(double [4][4], int);
 long double lg_pow(long double, int);



More information about the CIG-COMMITS mailing list