diff --git a/plugins/milk-extra-src/linalgebra/CMakeLists.txt b/plugins/milk-extra-src/linalgebra/CMakeLists.txt
index e4aaeb49..bd9cfd45 100644
--- a/plugins/milk-extra-src/linalgebra/CMakeLists.txt
+++ b/plugins/milk-extra-src/linalgebra/CMakeLists.txt
@@ -56,7 +56,9 @@ set(INCLUDEFILES
 	SingularValueDecomp_mkU.h
 )
 
-
+set(SCRIPTS
+    scripts/milk-modes-project
+)
 
 
 
@@ -145,3 +147,4 @@ message("---------------------------------------")
 
 install(TARGETS ${LIBNAME} DESTINATION lib)
 install(FILES ${INCLUDEFILES} DESTINATION include/${SRCNAME})
+install(PROGRAMS ${SCRIPTS} DESTINATION bin)
\ No newline at end of file
diff --git a/plugins/milk-extra-src/linalgebra/SGEMM.c b/plugins/milk-extra-src/linalgebra/SGEMM.c
index e69c6c29..c60bcf5e 100644
--- a/plugins/milk-extra-src/linalgebra/SGEMM.c
+++ b/plugins/milk-extra-src/linalgebra/SGEMM.c
@@ -76,7 +76,7 @@ static CLICMDARGDEF farg[] =
     {
         CLIARG_IMG,
         ".matB",
-        "output matrix B",
+        "input matrix B",
         "matA",
         CLIARG_VISIBLE_DEFAULT,
         (void **) &inmatB,
@@ -132,6 +132,7 @@ static CLICMDDATA CLIcmddata =
 static errno_t help_function()
 {
     printf("CPU or GPU matrix-matrix multiplication\n");
+    printf("Single precision only\n");
 
     return RETURN_SUCCESS;
 }
diff --git a/plugins/milk-extra-src/linalgebra/scripts/milk-modes-project b/plugins/milk-extra-src/linalgebra/scripts/milk-modes-project
new file mode 100755
index 00000000..bc3e4935
--- /dev/null
+++ b/plugins/milk-extra-src/linalgebra/scripts/milk-modes-project
@@ -0,0 +1,66 @@
+#!/usr/bin/env bash
+
+MSdescr="Linear projection of modesA onto space defined by modesB"
+
+MSextdescr="Project modesA (image cube) onto modesB (image cube)
+Neither modesA nor modesB are assumed to be orthogonal
+
+Main steps:
+- Perform PCA on modesB to build an orthonormal basis B
+- Decompose modesA according to basis B (decomposition coefficients are mcoeff)
+- Use mcoeff and basis B to reconstruct projection of modesA onto modesB
+- Compute the residual of the reconstruction = modesA - reconstruction
+
+Output:
+- projection
+- reconstruction
+"
+
+source milk-script-std-config
+
+RequiredCommands=( milk-all )
+RequiredFiles=()
+RequiredDirs=()
+
+
+MSarg+=( "inmodesA:string:input modes to be decomposed" )
+MSarg+=( "inmodesB:string:modes onto which to decompose" )
+MSarg+=( "SVDlim:float:SVD limit" )
+MSarg+=( "outimrec:string:reconstruction (projection)" )
+MSarg+=( "outimres:string:residual" )
+
+GPUindex="-1"
+MSopt+=( "g:gpu:setgpu:GPU[int]:GPU index, -1 for CPU" )
+function setgpu() {
+	GPUindex="$1"
+}
+
+
+source milk-argparse
+
+
+immA="${inputMSargARRAY[0]}"
+immB="${inputMSargARRAY[1]}"
+svdlim="${inputMSargARRAY[2]}"
+outimrec="${inputMSargARRAY[3]}"
+outimres="${inputMSargARRAY[4]}"
+
+milk-all << EOF
+loadfits "$immA" modesA
+loadfits "$immB" modesB
+linalg.compSVD .GPUdevice ${GPUindex}
+linalg.compSVD modesB svdU svdS svdV ${svdlim}
+listim
+linalg.sgemm .GPUdevice ${GPUindex}
+linalg.sgemm .transpA 1
+linalg.sgemm modesA svdU mcoeff
+linalg.sgemm .transpA 0
+linalg.sgemm .transpB 1
+linalg.sgemm svdU mcoeff imrec
+saveFITS imrec "$outimrec"
+imres=modesA-imrec
+saveFITS imres "$outimres"
+exitCLI
+EOF
+
+
diff --git a/scripts/milk-images-merge b/scripts/milk-images-merge
index 932472c9..b9bbaf86 100755
--- a/scripts/milk-images-merge
+++ b/scripts/milk-images-merge
@@ -53,7 +53,7 @@ while IFS= read -r filename; do
       echo "loadfits \"${filename}\" inimmerge${i}" >> $SF
       echo "immerge mergedim inimmerge${i} mergedim1 ${mergeaxis}" >> $SF
       echo "rm mergedim" >> $SF
-      echo "mv mergedim1 mergedim" >> $SF      
+      echo "mv mergedim1 mergedim" >> $SF
     fi
     ((i=i+1))
     initmerge=1