1 files changed, 74 insertions, 0 deletions

diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
index 310f6ee..0dccea2 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
@@ -814,3 +814,77 @@ int lu_l_C(KCMAT(a), DVEC(ipiv), CMAT(r)) {
     free(auxipiv);
     OK
 }
+
+////////////////////////////////////////////////////////////
+
+int multiplyR(KDMAT(a),KDMAT(b),DMAT(r)) {
+    REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
+    int i,j,k;
+    for (i=0;i<ar;i++) {
+        for(j=0;j<bc;j++) {
+            double temp = 0;
+            for(k=0;k<ac;k++) {
+                temp += ap[i*ac+k]*bp[k*bc+j];
+            }
+            rp[i*rc+j] = temp;
+        }
+    }
+    OK
+}
+
+int multiplyC(KCMAT(a),KCMAT(b),CMAT(r)) {
+    REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
+    int i,j,k;
+    for (i=0;i<ar;i++) {
+        for(j=0;j<bc;j++) {
+            doublecomplex temp = {0,0};
+            for(k=0;k<ac;k++) {
+                doublecomplex aik = ((doublecomplex*)ap)[i*ac+k];
+                doublecomplex bkj = ((doublecomplex*)bp)[k*bc+j];
+                //double w = aik.r+aik.i+bkj.r+bkj.i;
+                //if (w>w) exit(1);
+                //printf("%d",w>w);
+                temp.r += aik.r * bkj.r - aik.i * bkj.i;
+                temp.i += aik.r * bkj.i + aik.i * bkj.r;
+                //printf("%f %f %f %f \n",aik.r,aik.i,bkj.r,bkj.i);
+                //printf("%f %f %f \n",w,temp.r,temp.i);
+
+            }
+            ((doublecomplex*)rp)[i*rc+j] = temp;
+            //printf("%f %f\n",temp.r,temp.i);
+        }
+    }
+    OK
+}
+
+void dgemm_(char *, char *, integer *, integer *, integer *,
+           double *, const double *, integer *, const double *,
+           integer *, double *, double *, integer *);
+
+void zgemm_(char *, char *, integer *, integer *, integer *,
+           doublecomplex *, const doublecomplex *, integer *, const doublecomplex *,
+           integer *, doublecomplex *, doublecomplex *, integer *);
+
+
+int multiplyR2(KDMAT(a),KDMAT(b),DMAT(r)) {
+    REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
+    double alpha = 1;
+    double beta = 0;
+    integer m = ar;
+    integer n = bc;
+    integer k = ac;
+    int i,j;
+    dgemm_("N","N",&m,&n,&k,&alpha,ap,&m,bp,&k,&beta,rp,&m);
+    OK
+}
+
+int multiplyC2(KCMAT(a),KCMAT(b),CMAT(r)) {
+    REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
+    integer m = ar;
+    integer n = bc;
+    integer k = ac;
+    doublecomplex alpha = {1,0};
+    doublecomplex beta = {0,0};
+    zgemm_("N","N",&m,&n,&k,&alpha,(doublecomplex*)ap,&m,(doublecomplex*)bp,&k,&beta,(doublecomplex*)rp,&m);
+    OK
+}