[cig-commits] r6816 - in cs/cigma/trunk: include/cigma src tests

[luis at geodynamics.org]

Author: luis
Date: 2007-05-09 16:54:14 -0700 (Wed, 09 May 2007)
New Revision: 6816

Added:
   cs/cigma/trunk/include/cigma/det.h
   cs/cigma/trunk/src/det.c
   cs/cigma/trunk/tests/test_det.c
Log:
Explicit formulas for 3x3 and 4x4 determinants


Added: cs/cigma/trunk/include/cigma/det.h
===================================================================
--- cs/cigma/trunk/include/cigma/det.h	2007-05-09 23:51:44 UTC (rev 6815)
+++ cs/cigma/trunk/include/cigma/det.h	2007-05-09 23:54:14 UTC (rev 6816)
@@ -0,0 +1,13 @@
+#ifndef __CIGMA_DET_H__
+#define __CIGMA_DET_H__
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+double det3x3(double *m);
+double det4x4(double *m);
+
+#ifdef __cplusplus
+}
+#endif
+#endif

Added: cs/cigma/trunk/src/det.c
===================================================================
--- cs/cigma/trunk/src/det.c	2007-05-09 23:51:44 UTC (rev 6815)
+++ cs/cigma/trunk/src/det.c	2007-05-09 23:54:14 UTC (rev 6816)
@@ -0,0 +1,106 @@
+#include <cigma/det.h>
+
+double det3x3(double *m)
+{
+    //
+    // |m[0] m[1] m[2]|
+    // |m[3] m[4] m[5]|
+    // |m[6] m[7] m[8]|
+    //
+    // = m[0]*|m[4] m[5]|
+    //        |m[7] m[8]|
+    //
+    // - m[1]*|m[3] m[5]|
+    //        |m[6] m[8]|
+    //
+    // + m[2]*|m[3] m[4]|
+    //        |m[6] m[7]|
+    //
+    //
+    return m[0]*(m[4]*m[8] - m[5]*m[7])
+         - m[1]*(m[3]*m[8] - m[5]*m[6])
+         + m[2]*(m[3]*m[7] - m[4]*m[6]);
+}
+
+double det4x4(double *m)
+{
+    // |m[ 0] m[ 1] m[ 2] m[ 3]|
+    // |m[ 4] m[ 5] m[ 6] m[ 7]|
+    // |m[ 8] m[ 9] m[10] m[11]|
+    // |m[12] m[13] m[14] m[15]|
+    //
+    //
+    // = m[0] * |m[ 5] m[ 6] m[ 7]|
+    //          |m[ 9] m[10] m[11]|
+    //          |m[13] m[14] m[15]|
+    //
+    // - m[1] * |m[ 4] m[ 6] m[ 7]|
+    //          |m[ 8] m[10] m[11]|
+    //          |m[12] m[14] m[15]|
+    //
+    // + m[2] * |m[ 4] m[ 5] m[ 7]|
+    //          |m[ 8] m[ 9] m[11]|
+    //          |m[12] m[13] m[15]|
+    //
+    // - m[3] * |m[ 4] m[ 5] m[ 6]|
+    //          |m[ 8] m[ 9] m[10]|
+    //          |m[12] m[13] m[14]|
+    //
+    //
+    // = m[0] * ( + m[5] * |m[10] m[11]|
+    //                     |m[14] m[15]|
+    //
+    //            - m[6] * |m[ 9] m[11]|
+    //                     |m[13] m[15]|
+    //
+    //            + m[7] * |m[ 9] m[10]|
+    //                     |m[13] m[14]| )
+    //
+    // - m[1] * ( + m[4] * |m[10] m[11]|
+    //                     |m[14] m[15]|
+    //
+    //            - m[6] * |m[ 8] m[11]|
+    //                     |m[12] m[15]|
+    //
+    //            + m[7] * |m[ 8] m[10]|
+    //                     |m[12] m[14]| )
+    //
+    // + m[2] * ( + m[4] * |m[ 9] m[11]|
+    //                     |m[13] m[15]|
+    //
+    //            - m[5] * |m[ 8] m[11]|
+    //                     |m[12] m[15]|
+    //
+    //            + m[7] * |m[ 8] m[ 9]|
+    //                     |m[12] m[13]| )
+    //
+    // - m[3] * ( + m[4] * |m[ 9] m[10]|
+    //                     |m[13] m[14]|
+    //
+    //            - m[5] * |m[ 8] m[10]|
+    //                     |m[12] m[14]|
+    //
+    //            + m[6] * |m[ 8] m[ 9]|
+    //                     |m[12] m[13]| )
+
+    double c[4];
+
+    c[0] = m[5]*(m[10]*m[15]-m[11]*m[14])
+         - m[6]*(m[ 9]*m[15]-m[11]*m[13])
+         + m[7]*(m[ 9]*m[14]-m[10]*m[13]);
+
+    c[1] = m[4]*(m[10]*m[15]-m[11]*m[14])
+         - m[6]*(m[ 8]*m[15]-m[11]*m[12])
+         + m[7]*(m[ 8]*m[14]-m[10]*m[12]);
+
+    c[2] = m[4]*(m[ 9]*m[15]-m[11]*m[13])
+         - m[5]*(m[ 8]*m[15]-m[11]*m[12])
+         + m[7]*(m[ 8]*m[13]-m[ 9]*m[12]);
+
+    c[3] = m[4]*(m[ 9]*m[14]-m[10]*m[13])
+         - m[5]*(m[ 8]*m[14]-m[10]*m[12])
+         + m[6]*(m[ 8]*m[13]-m[ 9]*m[12]);
+
+    return m[0]*c[0] - m[1]*c[1] + m[2]*c[2] - m[3]*c[3];
+}
+

Added: cs/cigma/trunk/tests/test_det.c
===================================================================
--- cs/cigma/trunk/tests/test_det.c	2007-05-09 23:51:44 UTC (rev 6815)
+++ cs/cigma/trunk/tests/test_det.c	2007-05-09 23:54:14 UTC (rev 6816)
@@ -0,0 +1,57 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <cigma/det.h>
+
+double relative(double a, double b)
+{
+    return 100 * (b - a)/a;
+}
+
+void display(char *lhs, double a, double b)
+{
+    printf("%s = %+g  (%g%%)\n", lhs, b-a, relative(a,b));
+}
+
+int main(int argc, char *argv[])
+{
+    // x = numpy.random.rand(3,3)
+    double x[3*3] = {
+        0.946938064293, 0.938422952181, 0.547150183861,
+        0.357024276576, 0.0915660318218, 0.744489960953,
+        0.254916722044, 0.401961600912, 0.899098289403
+    };
+
+    // y = numpy.random.rand(4,4)
+    double y[4*4] = {
+        0.749255681616, 0.601726829068, 0.322624932027, 0.552755245286,
+        0.512189205379, 0.565678439509, 0.36170329247, 0.648798303295,
+        0.751073452252, 0.234564764206, 0.336464581301, 0.81318932981,
+        0.887471111472, 0.346416322698, 0.169225404909, 0.436142450865
+    };
+
+    // z = 1000 * numpy.random.rand(4,4)
+    double z[4*4] = {
+        929.692592005, 373.574551898, 886.825332236, 762.855520112,
+        207.72419653, 912.45502321, 140.297093513, 411.368439554,
+        813.530640335, 127.759885901, 788.306182949, 531.528945614,
+        694.198635678, 321.436908736, 680.495936175, 132.123707609
+    };
+
+    // numpy.linalg.det
+    double x_det = -0.262805708233;
+    double y_det = 0.0043930480585;
+    double z_det = 1839051314.33;
+    
+
+    double x_det3x3, y_det4x4, z_det4x4;
+
+    x_det3x3 = det3x3(x);
+    y_det4x4 = det4x4(y);
+    z_det4x4 = det4x4(z);
+
+    display("det3x3(x) - det(x)", x_det, x_det3x3);
+    display("det4x4(y) - det(y)", y_det, y_det4x4);
+    display("det4x4(z) - det(z)", z_det, z_det4x4);
+
+    return EXIT_SUCCESS;
+}