summaryrefslogtreecommitdiff
path: root/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
diff options
context:
space:
mode:
Diffstat (limited to 'lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c')
-rw-r--r--lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c74
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)) {
814 free(auxipiv); 814 free(auxipiv);
815 OK 815 OK
816} 816}
817
818////////////////////////////////////////////////////////////
819
820int multiplyR(KDMAT(a),KDMAT(b),DMAT(r)) {
821 REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
822 int i,j,k;
823 for (i=0;i<ar;i++) {
824 for(j=0;j<bc;j++) {
825 double temp = 0;
826 for(k=0;k<ac;k++) {
827 temp += ap[i*ac+k]*bp[k*bc+j];
828 }
829 rp[i*rc+j] = temp;
830 }
831 }
832 OK
833}
834
835int multiplyC(KCMAT(a),KCMAT(b),CMAT(r)) {
836 REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
837 int i,j,k;
838 for (i=0;i<ar;i++) {
839 for(j=0;j<bc;j++) {
840 doublecomplex temp = {0,0};
841 for(k=0;k<ac;k++) {
842 doublecomplex aik = ((doublecomplex*)ap)[i*ac+k];
843 doublecomplex bkj = ((doublecomplex*)bp)[k*bc+j];
844 //double w = aik.r+aik.i+bkj.r+bkj.i;
845 //if (w>w) exit(1);
846 //printf("%d",w>w);
847 temp.r += aik.r * bkj.r - aik.i * bkj.i;
848 temp.i += aik.r * bkj.i + aik.i * bkj.r;
849 //printf("%f %f %f %f \n",aik.r,aik.i,bkj.r,bkj.i);
850 //printf("%f %f %f \n",w,temp.r,temp.i);
851
852 }
853 ((doublecomplex*)rp)[i*rc+j] = temp;
854 //printf("%f %f\n",temp.r,temp.i);
855 }
856 }
857 OK
858}
859
860void dgemm_(char *, char *, integer *, integer *, integer *,
861 double *, const double *, integer *, const double *,
862 integer *, double *, double *, integer *);
863
864void zgemm_(char *, char *, integer *, integer *, integer *,
865 doublecomplex *, const doublecomplex *, integer *, const doublecomplex *,
866 integer *, doublecomplex *, doublecomplex *, integer *);
867
868
869int multiplyR2(KDMAT(a),KDMAT(b),DMAT(r)) {
870 REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
871 double alpha = 1;
872 double beta = 0;
873 integer m = ar;
874 integer n = bc;
875 integer k = ac;
876 int i,j;
877 dgemm_("N","N",&m,&n,&k,&alpha,ap,&m,bp,&k,&beta,rp,&m);
878 OK
879}
880
881int multiplyC2(KCMAT(a),KCMAT(b),CMAT(r)) {
882 REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
883 integer m = ar;
884 integer n = bc;
885 integer k = ac;
886 doublecomplex alpha = {1,0};
887 doublecomplex beta = {0,0};
888 zgemm_("N","N",&m,&n,&k,&alpha,(doublecomplex*)ap,&m,(doublecomplex*)bp,&k,&beta,(doublecomplex*)rp,&m);
889 OK
890}