summaryrefslogtreecommitdiff
path: root/lib/Data/Packed/Internal/auxi.c
diff options
context:
space:
mode:
Diffstat (limited to 'lib/Data/Packed/Internal/auxi.c')
-rw-r--r--lib/Data/Packed/Internal/auxi.c90
1 files changed, 12 insertions, 78 deletions
diff --git a/lib/Data/Packed/Internal/auxi.c b/lib/Data/Packed/Internal/auxi.c
index 7f83bcf..04dc7ad 100644
--- a/lib/Data/Packed/Internal/auxi.c
+++ b/lib/Data/Packed/Internal/auxi.c
@@ -4,14 +4,9 @@
4#include <gsl/gsl_matrix.h> 4#include <gsl/gsl_matrix.h>
5#include <gsl/gsl_math.h> 5#include <gsl/gsl_math.h>
6#include <gsl/gsl_errno.h> 6#include <gsl/gsl_errno.h>
7#include <gsl/gsl_fft_complex.h>
8#include <gsl/gsl_eigen.h>
9#include <gsl/gsl_integration.h>
10#include <gsl/gsl_deriv.h>
11#include <gsl/gsl_poly.h>
12#include <gsl/gsl_multimin.h>
13#include <gsl/gsl_complex.h> 7#include <gsl/gsl_complex.h>
14#include <gsl/gsl_complex_math.h> 8#include <gsl/gsl_complex_math.h>
9#include <gsl/gsl_cblas.h>
15#include <string.h> 10#include <string.h>
16#include <stdio.h> 11#include <stdio.h>
17 12
@@ -118,78 +113,6 @@ int constantC(gsl_complex* pval, CVEC(r)) {
118} 113}
119 114
120 115
121int multiplyR(int ta, KRMAT(a), int tb, KRMAT(b),RMAT(r)) {
122 //printf("%d %d %d %d %d %d\n",ar,ac,br,bc,rr,rc);
123 //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
124 DEBUGMSG("multiplyR (gsl_blas_dgemm)");
125 KDMVIEW(a);
126 KDMVIEW(b);
127 DMVIEW(r);
128 int k;
129 for(k=0;k<rr*rc;k++) rp[k]=0;
130 int debug = 0;
131 if(debug) {
132 printf("---------------------------\n");
133 printf("%p: ",ap); for(k=0;k<ar*ac;k++) printf("%f ",ap[k]); printf("\n");
134 printf("%p: ",bp); for(k=0;k<br*bc;k++) printf("%f ",bp[k]); printf("\n");
135 printf("%p: ",rp); for(k=0;k<rr*rc;k++) printf("%f ",rp[k]); printf("\n");
136 }
137 int res = gsl_blas_dgemm(
138 ta?CblasTrans:CblasNoTrans,
139 tb?CblasTrans:CblasNoTrans,
140 1.0, M(a), M(b),
141 0.0, M(r));
142 if(debug) {
143 printf("--------------\n");
144 printf("%p: ",ap); for(k=0;k<ar*ac;k++) printf("%f ",ap[k]); printf("\n");
145 printf("%p: ",bp); for(k=0;k<br*bc;k++) printf("%f ",bp[k]); printf("\n");
146 printf("%p: ",rp); for(k=0;k<rr*rc;k++) printf("%f ",rp[k]); printf("\n");
147 }
148 CHECK(res,res);
149 OK
150}
151
152int multiplyC(int ta, KCMAT(a), int tb, KCMAT(b),CMAT(r)) {
153 //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
154 DEBUGMSG("multiplyC (gsl_blas_zgemm)");
155 KCMVIEW(a);
156 KCMVIEW(b);
157 CMVIEW(r);
158 int k;
159 gsl_complex alpha, beta;
160 GSL_SET_COMPLEX(&alpha,1.,0.);
161 GSL_SET_COMPLEX(&beta,0.,0.);
162 //double *TEMP = (double*)malloc(rr*rc*2*sizeof(double));
163 //gsl_matrix_complex_view T = gsl_matrix_complex_view_array(TEMP,rr,rc);
164 for(k=0;k<rr*rc;k++) rp[k]=beta;
165 //for(k=0;k<2*rr*rc;k++) TEMP[k]=0;
166 int debug = 0;
167 if(debug) {
168 printf("---------------------------\n");
169 printf("%p: ",ap); for(k=0;k<2*ar*ac;k++) printf("%f ",((double*)ap)[k]); printf("\n");
170 printf("%p: ",bp); for(k=0;k<2*br*bc;k++) printf("%f ",((double*)bp)[k]); printf("\n");
171 printf("%p: ",rp); for(k=0;k<2*rr*rc;k++) printf("%f ",((double*)rp)[k]); printf("\n");
172 //printf("%p: ",T); for(k=0;k<2*rr*rc;k++) printf("%f ",TEMP[k]); printf("\n");
173 }
174 int res = gsl_blas_zgemm(
175 ta?CblasTrans:CblasNoTrans,
176 tb?CblasTrans:CblasNoTrans,
177 alpha, M(a), M(b),
178 beta, M(r));
179 //&T.matrix);
180 //memcpy(rp,TEMP,2*rr*rc*sizeof(double));
181 if(debug) {
182 printf("--------------\n");
183 printf("%p: ",ap); for(k=0;k<2*ar*ac;k++) printf("%f ",((double*)ap)[k]); printf("\n");
184 printf("%p: ",bp); for(k=0;k<2*br*bc;k++) printf("%f ",((double*)bp)[k]); printf("\n");
185 printf("%p: ",rp); for(k=0;k<2*rr*rc;k++) printf("%f ",((double*)rp)[k]); printf("\n");
186 //printf("%p: ",T); for(k=0;k<2*rr*rc;k++) printf("%f ",TEMP[k]); printf("\n");
187 }
188 CHECK(res,res);
189 OK
190}
191
192
193int diagR(KRVEC(d),RMAT(r)) { 116int diagR(KRVEC(d),RMAT(r)) {
194 REQUIRES(dn==rr && rr==rc,BAD_SIZE); 117 REQUIRES(dn==rr && rr==rc,BAD_SIZE);
195 DEBUGMSG("diagR"); 118 DEBUGMSG("diagR");
@@ -215,3 +138,14 @@ int diagC(KCVEC(d),CMAT(r)) {
215 } 138 }
216 OK 139 OK
217} 140}
141
142int conjugate(KCVEC(x),CVEC(t)) {
143 REQUIRES(xn==tn,BAD_SIZE);
144 DEBUGMSG("conjugate");
145 int k;
146 for (k=0; k<xn; k++) {
147 tp[k].dat[0] = xp[k].dat[0];
148 tp[k].dat[1] = - xp[k].dat[1];
149 }
150 OK
151}