linear algebra for modal ordering

plugins/milk-extra-src/linARfilterPred/build_linPF.c

@@ -710,6 +710,7 @@ static errno_t compute_function()
             imgU,
             imgeval,
             imgevec,
+            0,
             (*SVDeps),
             10000000,
             GPUdev,

plugins/milk-extra-src/linalgebra/CMakeLists.txt

@@ -15,6 +15,7 @@ set(SOURCEFILES
 	GPU_loop_MultMat_setup.c
 	GPU_SVD_computeControlMatrix.c
 	GPUloadCmat.c
+	GramSchmidt.c
 	magma_compute_SVDpseudoInverse_SVD.c
 	magma_compute_SVDpseudoInverse.c
 	magma_MatMatMult_testPseudoInverse.c
@@ -43,6 +44,7 @@ set(INCLUDEFILES
 	GPU_loop_MultMat_setup.h
 	GPU_SVD_computeControlMatrix.h
 	GPUloadCmat.h
+	GramSchmidt.h
 	magma_compute_SVDpseudoInverse_SVD.h
 	magma_compute_SVDpseudoInverse.h
 	magma_MatMatMult_testPseudoInverse.h

plugins/milk-extra-src/linalgebra/GramSchmidt.c

@@ -0,0 +1,204 @@
+#include <math.h>
+
+#include "CommandLineInterface/CLIcore.h"
+#include "COREMOD_iofits/COREMOD_iofits.h"
+
+#include "COREMOD_tools/COREMOD_tools.h"
+
+#include "SGEMM.h"
+
+
+static char *inmodes;
+static long  fpi_inmodes;
+
+static char *outmodes;
+static long  fpi_outmodes;
+
+
+static int32_t *GPUdevice;
+static long     fpi_GPUdevice;
+
+
+
+static CLICMDARGDEF farg[] =
+{
+    {
+        CLIARG_IMG,
+        ".inmodes",
+        "input modes",
+        "inm",
+        CLIARG_VISIBLE_DEFAULT,
+        (void **) &inmodes,
+        &fpi_inmodes
+    },
+    {
+        CLIARG_STR,
+        ".outmodes",
+        "output modes",
+        "outm",
+        CLIARG_VISIBLE_DEFAULT,
+        (void **) &outmodes,
+        &fpi_outmodes
+    },
+    {
+        // using GPU (99 : no GPU, otherwise GPU device)
+        CLIARG_INT32,
+        ".GPUdevice",
+        "GPU device, 99 for CPU",
+        "-1",
+        CLIARG_HIDDEN_DEFAULT,
+        (void **) &GPUdevice,
+        &fpi_GPUdevice
+    }
+};
+
+
+static CLICMDDATA CLIcmddata =
+{
+    "GramSchmidt", "Gram-Schmidt process", CLICMD_FIELDS_DEFAULTS
+};
+
+// detailed help
+static errno_t help_function()
+{
+    printf("Run Gram-Schmodt process\n");
+
+    return RETURN_SUCCESS;
+}
+
+
+errno_t GramSchmidt(
+    IMGID imginm,
+    IMGID *imgoutm,
+    int GPUdev
+)
+{
+    DEBUG_TRACE_FSTART();
+
+
+    // Compute cross product on input
+    IMGID imginxp  = mkIMGID_from_name("_outxp");
+    computeSGEMM(imginm, imginm, &imginxp, 1, 0, GPUdev);
+
+
+    // Create output
+    //
+    imcreatelikewiseIMGID(
+        imgoutm,
+        &imginm
+    );
+
+    uint32_t zsize;
+    uint32_t xysize = imginm.md->size[0];
+    if(imginm.md->naxis == 3)
+    {
+        zsize = imginm.md->size[2];
+        xysize *= imginm.md->size[1];
+    }
+    else
+    {
+        zsize = imginm.md->size[1];
+    }
+
+    printf("xysize = %u, zsize = %u\n", xysize, zsize);
+
+    for ( uint32_t kk=0; kk<zsize; kk++ )
+    {
+        // initializatoin
+        memcpy( &imgoutm->im->array.F[kk*xysize], &imginm.im->array.F[kk*xysize], sizeof(float)*xysize);
+
+        for ( uint32_t kk0 = 0; kk0 < kk; kk0++ )
+        {
+            // cross-product
+            double xpval = 0.0;
+
+            // square sum v0
+            double sqrsum0 = 0.0;
+
+            // square sum v1
+            double sqrsum1= 0.0;
+
+            for( uint32_t ii=0; ii<xysize; ii++)
+            {
+                float v0 = imgoutm->im->array.F[ kk0*xysize + ii];
+                float v1 = imgoutm->im->array.F[ kk*xysize + ii];
+
+                xpval += v0*v1;
+                sqrsum0 += v0*v0;
+                sqrsum1 += v1*v1;
+            }
+
+            float vcoeff = xpval / sqrsum0;
+
+            printf("  %5u  %5u   %f\n", kk, kk0, vcoeff);
+
+            for( uint32_t ii=0; ii<xysize; ii++)
+            {
+                imgoutm->im->array.F[ kk*xysize + ii] -= vcoeff * imgoutm->im->array.F[ kk0*xysize + ii];
+            }
+        }
+
+    }
+
+
+    list_image_ID();
+
+    DEBUG_TRACE_FEXIT();
+    return RETURN_SUCCESS;
+}
+
+
+
+
+
+
+static errno_t compute_function()
+{
+    DEBUG_TRACE_FSTART();
+
+    IMGID imginm = mkIMGID_from_name(inmodes);
+    resolveIMGID(&imginm, ERRMODE_ABORT);
+
+
+    IMGID imgoutm  = mkIMGID_from_name(outmodes);
+
+
+    INSERT_STD_PROCINFO_COMPUTEFUNC_INIT
+
+
+    INSERT_STD_PROCINFO_COMPUTEFUNC_LOOPSTART
+    {
+
+
+        GramSchmidt(imginm, &imgoutm, *GPUdevice);
+
+
+    }
+    INSERT_STD_PROCINFO_COMPUTEFUNC_END
+
+
+
+    DEBUG_TRACE_FEXIT();
+    return RETURN_SUCCESS;
+}
+
+
+
+
+INSERT_STD_FPSCLIfunctions
+
+
+
+
+// Register function in CLI
+errno_t
+CLIADDCMD_linalgebra__GramSchmidt()
+{
+
+    //CLIcmddata.FPS_customCONFsetup = customCONFsetup;
+    //CLIcmddata.FPS_customCONFcheck = customCONFcheck;
+    INSERT_STD_CLIREGISTERFUNC
+
+    return RETURN_SUCCESS;
+}
+

plugins/milk-extra-src/linalgebra/GramSchmidt.h

@@ -0,0 +1,13 @@
+#ifndef LINALGEBRA_GRAMSCHMIDT_H
+#define LINALGEBRA_GRAMSCHMIDT_H
+
+
+errno_t GramSchmidt(
+    IMGID imginm,
+    IMGID *imgoutm,
+    int GPUdev
+);
+
+errno_t CLIADDCMD_linalgebra__GramSchmidt();
+
+#endif

plugins/milk-extra-src/linalgebra/SGEMM.c

@@ -253,7 +253,15 @@ errno_t computeSGEMM(
         Ndim1 = inB_Mdim1;
     }
 
-    //printf("T %d %d  -> SGEMM  M=%d, N=%d, K=%d\n", TranspA, TranspB, Mdim, Ndim, Kdim);
+    printf("T %d %d  -> SGEMM  M=%d,(%d %d)  N=%d, (%d %d) K=%d\n",
+           TranspA, TranspB, Mdim, Mdim0, Mdim1, Ndim, Ndim0, Ndim1, Kdim);
+    printf("INPUT A  M %5d (%5d %5d)   x N %5d (%5d %5d)\n",
+           inA_Mdim, inA_Mdim0, inA_Mdim1,
+           inA_Ndim, inA_Ndim0, inA_Ndim1);
+    printf("INPUT B  M %5d (%5d %5d)   x N %5d (%5d %5d)\n",
+           inB_Mdim, inB_Mdim0, inB_Mdim1,
+           inB_Ndim, inB_Ndim0, inB_Ndim1);
+
 
 
     // Create output
@@ -277,6 +285,9 @@ errno_t computeSGEMM(
         outimg->size[1] = Mdim1;
         outimg->size[2] = outNdim;
     }
+
+    printf("OUTPUT  M %d   N %d  (%d %d %d)\n", outMdim, outNdim, outimg->size[0], outimg->size[1], outimg->size[2]);
+
     outimg->datatype = _DATATYPE_FLOAT;
     createimagefromIMGID(outimg);
 

plugins/milk-extra-src/linalgebra/SingularValueDecomp.c

@@ -47,6 +47,12 @@ static long  fpi_outS;
 static char *outV;
 static long  fpi_outV;
 
+// if V is 3D, set Vdim0 to its size[0]
+// otherwise leave at 0
+static uint32_t *Vdim0;
+static long   fpi_Vdim0;
+
+
 static float *svdlim;
 static long   fpi_svdlim;
 
@@ -103,6 +109,15 @@ static CLICMDARGDEF farg[] =
         (void **) &outV,
         &fpi_outV
     },
+    {
+        CLIARG_UINT32,
+        ".Vdim0",
+        "first dimension of V if 3D, 0 if 2D",
+        "0",
+        CLIARG_HIDDEN_DEFAULT,
+        (void **) &Vdim0,
+        &fpi_Vdim0
+    },
     {
         // Singular Value Decomposition limit
         CLIARG_FLOAT32,
@@ -230,6 +245,7 @@ errno_t compute_SVD(
     IMGID    imgU,
     IMGID    imgS,
     IMGID    imgV,
+    uint32_t Vdim0,
     float    SVlimit,
     uint32_t SVDmaxNBmode,
     int      GPUdev,
@@ -421,9 +437,19 @@ errno_t compute_SVD(
 
             if( imgV.ID == -1)
             {
-                imgV.naxis = 2;
-                imgV.size[0] = inNdim;
-                imgV.size[1] = NBmode; //inNdim;
+                if( Vdim0 == 0)
+                {
+                    imgV.naxis = 2;
+                    imgV.size[0] = inNdim;
+                    imgV.size[1] = NBmode; //inNdim;
+                }
+                else
+                {
+                    imgV.naxis = 3;
+                    imgV.size[0] = Vdim0;
+                    imgV.size[1] = inNdim/Vdim0;
+                    imgV.size[2] = NBmode; //inNdim;
+                }
                 createimagefromIMGID(&imgV);
             }
         }
@@ -677,7 +703,7 @@ static errno_t compute_function()
     {
 
 
-        compute_SVD(imginM, imgU, imgS, imgV, *svdlim, *maxNBmode, *GPUdevice, *compmode);
+        compute_SVD(imginM, imgU, imgS, imgV, *Vdim0, *svdlim, *maxNBmode, *GPUdevice, *compmode);
 
 
     }

plugins/milk-extra-src/linalgebra/SingularValueDecomp.h

@@ -13,6 +13,7 @@ errno_t compute_SVD(
     IMGID imgU,
     IMGID imgeigenval,
     IMGID imgV,
+    uint32_t Vdim0,
     float SVlimit,
     uint32_t SVDmaxNBmode,
     int GPUdev,