1 files changed, 405 insertions, 0 deletions

diff --git a/SD-VBS/common/toolbox/MultiNcut/a_times_b_cmplx.c b/SD-VBS/common/toolbox/MultiNcut/a_times_b_cmplx.c
new file mode 100755
index 0000000..25def92
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/a_times_b_cmplx.c
@@ -0,0 +1,405 @@
+/*================================================================
+a_times_b_cmplx.c = used by a couple of mex functions 
+provide Matrix vector multiplications,
+and solve triangular systems
+(sparse matrix and full vector)
+
+CSC_CmplxVecMult_CAB_double, CSR_CmplxVecMult_CAB_double, 
+CSCsymm_CmplxVecMult_CAB_double added by Mirko Visontai (10/24/2003)
+
+*=================================================================*/
+# include "math.h"
+
+
+/*C<-a*A*B+C*/
+void CSC_VecMult_CaABC_double(
+                 const int m,  const int k, const double alpha,
+                 const double *val, const int *indx, 
+                 const int *pntrb,
+                 const double *b,  
+                 double *c)
+{
+  int i,j,jb,je;
+
+    for (i=0;i!=k;i++){
+      jb = pntrb[i];
+      je = pntrb[i+1];
+      for (j=jb;j!=je;j++)
+        c[indx[j]] += alpha * b[i] * val[j];
+    }
+}
+
+/*C<-a*A'*B+C*/
+void CSR_VecMult_CaABC_double(
+                 const int k,  const int m, const double alpha,
+                 const double *val, const int *indx, 
+                 const int *pntrb,
+                 const double *b,  
+                 double *c)
+{
+  double t;
+  const double *pval;
+  int i,j,jb,je;
+
+    pval = val;
+    for (i=0;i!=m;i++) {
+      t = 0;
+      jb = pntrb[i];
+      je = pntrb[i+1];
+      for (j=jb;j!=je;j++)
+        t += alpha * b[indx[j]] * (*pval++);
+      c[i] += t;
+    }
+}
+
+
+/*C<-A*b */
+void CSC_VecMult_CAB_double(
+                 const int m,  const int k, /*nb_rows, nb_columns*/
+                 const double *val, const int *indx, 
+                 const int *pntrb,                  
+                 const double *b,  
+                 double *c                 
+                 )
+{
+  int i,j,jb,je;
+  double *pc=c;                                       
+    for (i=0;i!=m;i++) *pc++ = 0;                 
+
+    for (i=0;i!=k;i++){
+      jb = pntrb[i];
+      je = pntrb[i+1];
+      for (j=jb;j!=je;j++)
+        c[indx[j]] +=  b[i] * val[j];
+    }
+}
+
+/*C<-A*b (complex)*/
+void CSC_CmplxVecMult_CAB_double(
+                 const int m,  const int k, 
+                 const double *valr, const double *vali,
+                                 const int *indx, 
+                 const int *pntrb,                  
+                 const double *br, const double *bi,  
+                 double *cr, double *ci                 
+                 )
+{
+  int i,j,jb,je;
+  double *pcr=cr;                                       
+  double *pci=ci;                                       
+    for (i=0;i!=m;i++){
+                *pcr++ = 0.0;
+                *pci++ = 0.0;
+        }
+
+    for (i=0;i!=k;i++){
+      jb = pntrb[i];
+      je = pntrb[i+1];
+      for (j=jb;j!=je;j++){
+        cr[indx[j]] +=  (br[i] * valr[j]) - (bi[i] * vali[j]);
+        ci[indx[j]] +=  (br[i] * vali[j]) + (bi[i] * valr[j]);
+          }
+    }
+}
+
+/*C<-A'*b 
+ plus rapide que CSC_VecMult_CAB_double */
+void CSR_VecMult_CAB_double(
+                 const int k,  const int m,
+                 const double *val, const int *indx, 
+                 const int *pntrb,
+                 const double *b,  
+                 double *c
+                 )
+{
+  double t;
+  const double *pval;
+  double *pc=c;                                       
+  int i,j,jb,je;
+ 
+    for (i=0;i!=m;i++) *pc++ = 0;                 
+
+    pval = val;
+    for (i=0;i!=m;i++) {
+      t = 0;
+      jb = pntrb[i];
+      je = pntrb[i+1];
+      for (j=jb;j!=je;j++)
+        t +=  b[indx[j]] * (*pval++);
+      c[i] += t;
+    }
+}
+
+/*C<-A'*b (complex) 
+ plus rapide que CSC_VecMult_CAB_double */
+void CSR_CmplxVecMult_CAB_double(
+                 const int k,  const int m, 
+                 const double *valr, const double *vali,
+                                 const int *indx, 
+                 const int *pntrb,
+                 const double *br, const double *bi,  
+                 double *cr, double *ci
+                 )
+{
+  double tr, ti;
+  const double *pvalr;
+  const double *pvali;
+  double *pcr=cr;                                       
+  double *pci=ci;                                       
+  int i,j,jb,je;
+ 
+    for (i=0;i!=m;i++){
+                *pcr++ = 0.0;                 
+                *pci++ = 0.0;
+        }               
+
+    pvalr = valr;
+    pvali = vali;
+    for (i=0;i!=m;i++) {
+      tr = 0.0;
+          ti = 0.0;
+      jb = pntrb[i];
+      je = pntrb[i+1];
+      for (j=jb;j!=je;j++){
+        tr +=  (br[indx[j]] * (*pvalr)) - (bi[indx[j]] * (*pvali));
+        ti +=  (br[indx[j]] * (*pvali++)) + (bi[indx[j]] * (*pvalr++));
+          }
+      cr[i] += tr;
+          ci[i] += ti;
+    }
+}
+
+
+
+/* C<-A*b (A is symmetric) */
+void CSRsymm_VecMult_CAB_double(
+                 const int k,  const int m, 
+                 const double *val, const int *indx, 
+                 const int *pntrb,  
+                 const double *b,  
+                 double *c 
+                 )
+{
+  const double *pval;
+  double *pc=c;                                       
+  int i,j;
+  int jj;
+  int rpntrb, rpntre;
+  int index, nvals;
+  
+      
+    for (i=0;i!=m;i++) *pc++ = 0;                 
+    pval = val;
+    for (j=0;j!=k;j++){
+      rpntrb = pntrb[j];
+      rpntre = pntrb[j+1];
+      for (jj=rpntrb;jj!=rpntre;jj++) {
+        index = indx[jj];
+        if ( index == j ) {
+          c[j] +=  b[j] * (*pval++);
+          continue;
+        }
+                if ( index > j ) {
+                c[index] +=  b[j] * (*pval);
+
+                c[j] +=  b[index] * (*pval++);
+                }
+                else {
+                        pval++;
+                }
+      }
+    }
+}
+
+
+/* C<-A*b (A is symmetric and complex) */
+void CSRsymm_CmplxVecMult_CAB_double(
+                 const int k,  const int m, 
+                 const double *valr, const double *vali,
+                                 const int *indx, 
+                 const int *pntrb,  
+                 const double *br, const double *bi,  
+                 double *cr, double *ci
+                 )
+{
+  const double *pvalr, *pvali;
+  double *pcr=cr;                                       
+  double *pci=ci;                                       
+  int i,j;
+  int jj;
+  int rpntrb, rpntre;
+  int index, nvals;
+  
+      
+    for (i=0;i!=m;i++){
+                *pcr++ = 0.0;                 
+                *pci++ = 0.0;                 
+        }
+
+    pvalr = valr;
+    pvali = vali;
+    for (j=0;j!=k;j++){
+      rpntrb = pntrb[j];
+      rpntre = pntrb[j+1];
+      for (jj=rpntrb;jj!=rpntre;jj++) {
+        index = indx[jj];
+        if ( index == j ) {
+          cr[j] +=  (br[j] * (*pvalr)) - (bi[j] * (*pvali));
+          ci[j] +=  (br[j] * (*pvali++)) + (bi[j] * (*pvalr++));
+          continue;
+        }
+        if ( index > j ) {
+                cr[index] +=  (br[j] * (*pvalr)) - (bi[j] * (*pvali));
+                ci[index] +=  (br[j] * (*pvali)) + (bi[j] * (*pvalr));
+        
+                cr[j] +=  (br[index] * (*pvalr)) - (bi[index] * (*pvali));
+                ci[j] +=  (br[index] * (*pvali++)) + (bi[index] * (*pvalr++));
+                }
+                else {
+                        pvalr++;
+                        pvali++;
+                }
+        
+      }
+    }
+}
+
+
+/*C<-A\B; with Lower triangular A*/
+void CSC_VecTriangSlvLD_CAB_double( 
+                 const int m,  
+                 const double *val,
+                 const int *indx, const int *pntrb, 
+                 const double *b,  
+                 double *c)
+{
+  int i, j, jb, je;
+  double *pc=c;
+  double z; 
+
+    for (i=0;i!=m;i++){
+      *pc = b[i];
+      pc++;
+    }                                     
+
+    pc=c;
+    for (i=0;i!=m;i++) {
+      jb = pntrb[i];
+      je = pntrb[i+1];
+      z =  pc[i] / val[jb];
+      pc[i] = z;
+      for (j=jb+1;j<je;j++) {
+        c[indx[j]] -= z*val[j];
+      }
+    }
+}
+
+/*C<-A\B; with Upper triangular A*/
+void CSC_VecTriangSlvUD_CAB_double(
+                 const int m,  
+                   const double *val,
+                 const int *indx, const int *pntrb,
+                 const double *b,  
+                 double *c)
+{
+  int i, j, jb, je, index;
+  double *pc=c;
+  double z; 
+
+    for (i=0;i!=m;i++){
+      *pc = b[i];
+      pc++;
+    }                                     
+
+    pc=c;
+    for (i=m-1;i!=-1;i--) {
+      jb = pntrb[i];
+      je = pntrb[i+1]-1;
+      z = pc[i] /val[je];
+      pc[i] = z;
+      for (j=jb;j<je;j++) {
+          c[indx[j]] -= z * val[j];
+      }
+    }
+}
+/*C<-A'\B; where A is upper (little slower than CSC)*/
+void CSR_VecTriangSlvLD_CAB_double(
+                 const int m,   
+                 const double *val, 
+                 const int *indx, const int *pntrb,
+                 const double *b,   
+                 double *c)
+{
+  int i, j, jb, je, index;
+  double *pc=c;
+  double z; 
+  double valtmp;
+
+    pc=c;
+    for (i=0;i!=m;i++) {
+      z = 0;
+      jb = pntrb[i];
+      je = pntrb[i+1];
+      for (j=jb;j<je;j++) {
+        index = indx[j];
+        if ( index == i ) {
+          valtmp = val[j];
+        } else {
+          z += c[index] * val[j];
+        }
+      }
+      pc[i] = (b[i] - z) / valtmp;
+    }
+}
+
+/*C<-A'\B; where A is lower  (little slower than CSC)*/
+void CSR_VecTriangSlvUD_CAB_double(
+                 const int m,   
+                   const double *val, 
+                 const int *indx, const int *pntrb, 
+                 const double *b,   
+                 double *c)
+{
+  int i, j, jb, je, index;
+  double *pc=c;
+  double valtmp;
+  double z; 
+
+    pc=c;
+    for (i=m-1;i!=-1; i--) {
+      z = 0;
+      jb = pntrb[i];
+      je =  pntrb[i+1];
+      for (j=jb+1; j<je; j++) {
+        z += c[indx[j]] * val[j];
+      }
+      pc[i] = (b[i] - z) / val[jb];
+    }
+}
+
+/*C<-A*B, where A is (m,k) and B is (k,n)*/
+void CSC_MatMult_CAB_double(
+                 const int m, const int n, const int k, 
+                 const double *val, const int *indx, 
+                 const int *pntrb,
+                 const double *b, const int ldb, 
+                 double *c, const int ldc)
+{
+  int i,j,jb,je;
+  double *pc=c;                                       
+  int l;                                               
+
+  for (l=0;l!=n;l++)                           
+    for (i=0;i!=m;i++) *pc++ = 0;                 
+
+  for (l=0;l!=n;l++) {                                 
+    for (i=0;i!=k;i++){
+      jb = pntrb[i];
+      je = pntrb[i+1];
+      for (j=jb;j!=je;j++)
+        c[indx[j]] +=  b[i] * val[j];
+    }
+    /*c += ldc; b += ldb;  */                              
+    c += m; b += m;                                
+  }                                                    
+}