add co-input

milk-org · Nov 14, 2023 · ac40c91 · ac40c91
1 parent 6aca545
commit ac40c91
Show file tree

Hide file tree

Showing 2 changed files with 23 additions and 6 deletions.
diff --git a/plugins/milk-extra-src/linalgebra/GramSchmidt.c b/plugins/milk-extra-src/linalgebra/GramSchmidt.c
@@ -95,8 +95,9 @@ errno_t GramSchmidt(
 
 
     // Compute cross product on input
-    IMGID imginxp  = mkIMGID_from_name("_outxp");
-    computeSGEMM(imginm, imginm, &imginxp, 1, 0, GPUdev);
+    //
+    //IMGID imginxp  = mkIMGID_from_name("_outxp");
+    //computeSGEMM(imginm, imginm, &imginxp, 1, 0, GPUdev);
 
 
     // Create output
@@ -118,11 +119,25 @@ errno_t GramSchmidt(
         zsize = imginm.md->size[1];
     }
 
+    uint32_t xysizeaux = 0;
+    if ( imgaux.ID != -1)
+    {
+        xysizeaux = imgaux.md->size[0];
+        if(imginm.md->naxis == 3)
+        {
+            xysizeaux *= imgaux.md->size[1];
+        }
+    }
+
+
     printf("xysize = %u, zsize = %u\n", xysize, zsize);
 
+    printf("\n");
     for ( uint32_t kk=0; kk<zsize; kk++ )
     {
-        // initializatoin
+        printf("\rGS mode %6u / %6u     ", kk, zsize);
+
+        // initialization
         memcpy( &imgoutm->im->array.F[kk*xysize], &imginm.im->array.F[kk*xysize], sizeof(float)*xysize);
 
         for ( uint32_t kk0 = 0; kk0 < kk; kk0++ )
@@ -157,13 +172,14 @@ errno_t GramSchmidt(
 
             if ( imgaux.ID != -1)
             {
-                for( uint32_t ii=0; ii<xysize; ii++)
+                for( uint32_t ii=0; ii<xysizeaux; ii++)
                 {
-                    imgaux.im->array.F[ kk*xysize + ii] -= vcoeff * imgaux.im->array.F[ kk0*xysize + ii];
+                    imgaux.im->array.F[ kk*xysizeaux + ii] -= vcoeff * imgaux.im->array.F[ kk0*xysizeaux + ii];
                 }
             }
         }
     }
+    printf("\n");
 
 
     list_image_ID();

diff --git a/plugins/milk-extra-src/linalgebra/SingularValueDecomp.c b/plugins/milk-extra-src/linalgebra/SingularValueDecomp.c
@@ -696,9 +696,10 @@ errno_t compute_SVD(
         for(int kk=0; kk<imgvnmodes.size[lastaxis]; kk++)
         {
             float mfact = imgS.im->array.F[kk];
+            //printf("mfact %4d = %f\n", kk, mfact);
             for(long ii=0; ii<framesize; ii++)
             {
-                imgvnmodes.im->array.F[kk*framesize+ii] = imgV.im->array.F[kk*framesize+ii] / mfact;
+                imgvnmodes.im->array.F[kk*framesize+ii] = imgV.im->array.F[kk*framesize+ii] * mfact;
             }
         }