00001 #include "operator.hpp"
00002 #include <math.h>
00003
00004 SmoluchowskiCalc::SmoluchowskiCalc(const Skeleton& K_skel,
00005 const Skeleton& Psi_skel, double _H, unsigned _N)
00006 : K(K_skel), Psi(Psi_skel), N(_N), H(_H), h(H / double(N))
00007 { }
00008
00009
00010
00011
00012 SmoluchowskiCalcDirect::SmoluchowskiCalcDirect(const Skeleton& K_skel,
00013 const Skeleton& Psi_skel, double _H, unsigned _N)
00014 : SmoluchowskiCalc(K_skel, Psi_skel, _H, _N), Uf(0)
00015 { }
00016
00017 SmoluchowskiCalcDirect::~SmoluchowskiCalcDirect()
00018 {
00019 if(Uf) {
00020 delete [] Uf;
00021 }
00022 }
00023
00024
00025 void SmoluchowskiCalcDirect::calc_L1(const double *f,
00026 const double *g, double *res)
00027 {
00028 unsigned i, j, k, r;
00029 r = K.GetRank();
00030 for(i = 0; i <= N; ++i) {
00031 res[i] = 0;
00032 for(k = 0; k < r; ++k) {
00033 for(j = 0; j <= i; ++j) {
00034 res[i] += K.GetU().elem(i - j, k) * f[i - j] *
00035 K.GetV().elem(j, k) * g[j];
00036 }
00037 res[i] -= 0.5 * (K.GetU().elem(i, k) * f[i] *
00038 K.GetV().elem(0, k) * g[0] +
00039 K.GetU().elem(0, k) * f[0] *
00040 K.GetV().elem(i, k) * g[i]);
00041 }
00042 res[i] *= h;
00043 }
00044 }
00045
00046 void SmoluchowskiCalcDirect::calc_L2(const double *f,
00047 const double *g, double *res)
00048 {
00049 double tmp;
00050 unsigned i, j, k, r;
00051 const Matrix& U = K.GetU();
00052 const Matrix& V = K.GetV();
00053 r = K.GetRank();
00054 for(i = 0; i <= N; ++i) {
00055 res[i] = 0;
00056 for(k = 0; k < r; ++k) {
00057 tmp = 0;
00058 for(j = i; j <= N; ++j) {
00059 tmp += U.elem(j - i, k) * f[j - i] * g[j];
00060 }
00061 tmp -= 0.5 * (U.elem(0, k) * f[0] * g[i] +
00062 U.elem(N - i, k) * f[N - i] * g[N]);
00063 res[i] += V.elem(i, k) * tmp;
00064 }
00065 res[i] *= h;
00066 }
00067 }
00068
00069 void SmoluchowskiCalcDirect::fix_f_arg(const double *f)
00070 {
00071 unsigned i, k, r;
00072 r = K.GetRank();
00073 if(!Uf) {
00074 Uf = new double[(N + 1) * r];
00075 }
00076 for(k = 0; k < r; ++k) {
00077 for(i = 0; i <= N; ++i) {
00078 Uf[(N + 1) * k + i] = K.GetU().elem(i, k) * f[i];
00079 }
00080 }
00081 }
00082
00083 void
00084 SmoluchowskiCalcDirect::calc_L1_fix_f(const double *g, double *res)
00085 {
00086 unsigned i, j, k, r;
00087 #if HANDLEXCEPTION
00088 if(!Uf) {
00089 throw "f must be set by fix_f_arg";
00090 }
00091 #endif
00092 r = K.GetRank();
00093 for(i = 0; i <= N; ++i) {
00094 res[i] = 0;
00095 for(k = 0; k < r; ++k) {
00096 for(j = 0; j <= i; ++j) {
00097 res[i] += Uf[(N + 1) * k + i - j] *
00098 K.GetV().elem(j, k) * g[j];
00099 }
00100 res[i] -= 0.5 * (Uf[(N + 1) * k + i] *
00101 K.GetV().elem(0, k) * g[0] +
00102 Uf[(N + 1) * k] *
00103 K.GetV().elem(i, k) * g[i]);
00104 }
00105 res[i] *= h;
00106 }
00107 }
00108
00109 void
00110 SmoluchowskiCalcDirect::calc_L2_fix_f(const double *g, double *res)
00111 {
00112 double tmp;
00113 unsigned i, j, k, r;
00114 #if HANDLEXCEPTION
00115 if(!Uf) {
00116 throw "f must be set by fix_f_arg";
00117 }
00118 #endif
00119 r = K.GetRank();
00120 for(i = 0; i <= N; ++i) {
00121 res[i] = 0;
00122 for(k = 0; k < r; ++k) {
00123 tmp = 0;
00124 for(j = i; j <= N; ++j) {
00125 tmp += Uf[(N + 1) * k + j - i] * g[j];
00126 }
00127 tmp -= 0.5 * (Uf[(N + 1) * k] * g[i] +
00128 Uf[(N + 1) * k + N - i] * g[N]);
00129 res[i] += K.GetV().elem(i, k) * tmp;
00130 }
00131 res[i] *= h;
00132 }
00133 }
00134
00135
00136
00137 void SmoluchowskiCalcDirect::calc_L3(const double *g, double *res) const
00138 {
00139 unsigned i, j;
00140 for(i = 0; i <= N; ++i) {
00141 res[i] = 0;
00142 for(j = 0; j <= i; ++j) {
00143 res[i] += Psi.GetM().elem(i, j) * g[j];
00144 }
00145 res[i] = h * (res[i] - 0.5 * (Psi.GetM().elem(i, 0) * g[0] +
00146 Psi.GetM().elem(i, i) * g[i]));
00147 }
00148 }
00149
00150 void SmoluchowskiCalcDirect::calc_L4(const double *g, double *res) const
00151 {
00152 unsigned i, j;
00153 for(i = 0; i <= N; ++i) {
00154 res[i] = 0;
00155 for(j = i; j <= N; ++j) {
00156 res[i] += Psi.GetM().elem(j, i) * g[j];
00157 }
00158 res[i] = h * (res[i] - 0.5 * (Psi.GetM().elem(i, i) * g[i] +
00159 Psi.GetM().elem(N, i) * g[N]));
00160 }
00161 }
00162
00163 void SmoluchowskiCalcDirect::calc_L5(const double *g, double *res) const
00164 {
00165 unsigned i, j;
00166 for(i = 0; i <= N; ++i) {
00167 res[i] = 0;
00168 for(j = 0; j <= N; ++j) {
00169 res[i] += K.GetM().elem(i, j) * g[j];
00170 }
00171 res[i] = h * (res[i] - 0.5 * (K.GetM().elem(i, 0) * g[0] +
00172 K.GetM().elem(i, N) * g[N]));
00173 }
00174 }
00175
00176
00177
00178
00179 SmoluchowskiCalcFast::SmoluchowskiCalcFast(const Skeleton& K_skel,
00180 const Skeleton& Psi_skel, double _H, unsigned _N)
00181 : SmoluchowskiCalc(K_skel, Psi_skel, _H, _N),
00182 N_log2(ceil(log2(2 * N + 1))), N_2(1 << N_log2), Fourier(N_log2),
00183 Uf(0), Uf_fourier(0)
00184 {
00185 unsigned r1, r2;
00186 r1 = K.GetRank();
00187 r2 = Psi.GetRank();
00188 tmp_mv = new double[r1 > r2 ? r1 : r2];
00189 tmp_complex0 = new Complex[N_2];
00190 tmp_complex1 = new Complex[N_2];
00191 tmp_complex2 = new Complex[N_2];
00192 }
00193
00194 SmoluchowskiCalcFast::~SmoluchowskiCalcFast()
00195 {
00196 delete [] tmp_mv;
00197 delete [] tmp_complex0;
00198 delete [] tmp_complex1;
00199 delete [] tmp_complex2;
00200 if(Uf) {
00201 delete [] Uf;
00202 delete [] Uf_fourier;
00203 }
00204 }
00205
00206
00207 void SmoluchowskiCalcFast::calc_L1(const double *f,
00208 const double *g, double *res)
00209 {
00210 unsigned i, k, r;
00211 r = K.GetRank();
00212 for(i = 0; i <= N; ++i) {
00213 res[i] = 0;
00214 }
00215 for(k = 0; k < r; ++k) {
00216 for(i = 0; i <= N; ++i) {
00217 tmp_complex0[i] = K.GetU().elem(i, k) * f[i];
00218 }
00219 Fourier.MatvecFast(tmp_complex0, tmp_complex1);
00220 for(i = 0; i <= N; ++i) {
00221 tmp_complex0[i] = K.GetV().elem(i, k) * g[i];
00222 }
00223 Fourier.MatvecFast(tmp_complex0, tmp_complex2);
00224 for(i = 0; i < N_2; ++i) {
00225 tmp_complex1[i] *= tmp_complex2[i];
00226 }
00227 Fourier.MatvecFast(tmp_complex1, tmp_complex2,
00228 FourierTransform::conj);
00229 for(i = 0; i <= N; ++i) {
00230 res[i] +=
00231 h * (tmp_complex2[i].GetRe() / double(N_2) -
00232 0.5 * (K.GetU().elem(i, k) * f[i] *
00233 tmp_complex0[0].GetRe() +
00234 K.GetU().elem(0, k) * f[0] *
00235 tmp_complex0[i].GetRe()));
00236 }
00237 }
00238 }
00239
00240 void SmoluchowskiCalcFast::calc_L2(const double *f,
00241 const double *g, double *res)
00242 {
00243 unsigned i, k, r;
00244 r = K.GetRank();
00245 for(i = 0; i <= N; ++i) {
00246 res[i] = 0;
00247 }
00248 for(k = 0; k < r; ++k) {
00249 for(i = 0; i <= N; ++i) {
00250 tmp_complex0[N_2 - 1 - i] = g[i];
00251 }
00252 Fourier.MatvecFast(tmp_complex0, tmp_complex1);
00253 for(i = 0; i <= N; ++i) {
00254 tmp_complex0[i] = K.GetU().elem(i, k) * f[i];
00255 }
00256 for(i = N + 1; i < N_2; ++i) {
00257 tmp_complex0[i] = 0;
00258 }
00259 Fourier.MatvecFast(tmp_complex0, tmp_complex2);
00260 for(i = 0; i < N_2; ++i) {
00261 tmp_complex1[i] *= tmp_complex2[i];
00262 }
00263 Fourier.MatvecFast(tmp_complex1, tmp_complex2,
00264 FourierTransform::conj);
00265 for(i = 0; i <= N; ++i) {
00266 res[i] += h * K.GetV().elem(i, k) *
00267 (tmp_complex2[N_2 - 1 - i].GetRe() / double(N_2) -
00268 0.5 * (tmp_complex0[0].GetRe() * g[i] +
00269 tmp_complex0[N - i].GetRe() * g[N]));
00270 }
00271 }
00272 }
00273
00274 void SmoluchowskiCalcFast::fix_f_arg(const double *f)
00275 {
00276 unsigned i, k, r;
00277 r = K.GetRank();
00278 if(!Uf) {
00279 Uf = new double[(N + 1) * K.GetRank()];
00280 Uf_fourier = new Complex[N_2 * r];
00281 }
00282 for(k = 0; k < r; ++k) {
00283 for(i = 0; i <= N; ++i) {
00284 Uf[(N + 1) * k + i] = K.GetU().elem(i, k) * f[i];
00285 tmp_complex1[i] = Uf[(N + 1) * k + i];
00286 }
00287 Fourier.MatvecFast(tmp_complex1, Uf_fourier + N_2 * k);
00288 }
00289 }
00290
00291 void
00292 SmoluchowskiCalcFast::calc_L1_fix_f(const double *g, double *res)
00293 {
00294 unsigned i, k, r;
00295 #if HANDLEXCEPTION
00296 if(!Uf) {
00297 throw "f must be set by fix_f_arg";
00298 }
00299 #endif
00300 r = K.GetRank();
00301 for(i = 0; i <= N; ++i) {
00302 res[i] = 0;
00303 }
00304 for(k = 0; k < r; ++k) {
00305 for(i = 0; i <= N; ++i) {
00306 tmp_complex0[i] = K.GetV().elem(i, k) * g[i];
00307 }
00308 Fourier.MatvecFast(tmp_complex0, tmp_complex1);
00309 for(i = 0; i < N_2; ++i) {
00310 tmp_complex1[i] *= Uf_fourier[N_2 * k + i];
00311 }
00312 Fourier.MatvecFast(tmp_complex1, tmp_complex2,
00313 FourierTransform::conj);
00314 for(i = 0; i <= N; ++i) {
00315 res[i] += h * (tmp_complex2[i].GetRe() / double(N_2) -
00316 0.5 * (Uf[(N + 1) * k + i] * tmp_complex0[0].GetRe() +
00317 Uf[(N + 1) * k] * tmp_complex0[i].GetRe()));
00318 }
00319 }
00320 }
00321
00322 void
00323 SmoluchowskiCalcFast::calc_L2_fix_f(const double *g, double *res)
00324 {
00325 unsigned i, k, r;
00326 #if HANDLEXCEPTION
00327 if(!Uf) {
00328 throw "f must be set by fix_f_arg";
00329 }
00330 #endif
00331 r = K.GetRank();
00332 for(i = 0; i <= N; ++i) {
00333 res[i] = 0;
00334 }
00335 for(k = 0; k < r; ++k) {
00336 for(i = 0; i <= N; ++i) {
00337 tmp_complex0[N_2 - 1 - i] = g[i];
00338 }
00339 Fourier.MatvecFast(tmp_complex0, tmp_complex1);
00340 for(i = 0; i < N_2; ++i) {
00341 tmp_complex1[i] *= Uf_fourier[N_2 * k + i];
00342 tmp_complex0[i] = 0;
00343 }
00344 Fourier.MatvecFast(tmp_complex1, tmp_complex2,
00345 FourierTransform::conj);
00346 for(i = 0; i <= N; ++i) {
00347 res[i] += h * K.GetV().elem(i, k) *
00348 (tmp_complex2[N_2 - 1 - i].GetRe() / double(N_2) -
00349 0.5 * (Uf[(N + 1) * k] * g[i] +
00350 Uf[(N + 1) * k + N - i] * g[N]));
00351 }
00352 }
00353 }
00354
00355
00356 void SmoluchowskiCalcFast::calc_L3(const double *g, double *res) const
00357 {
00358 unsigned i;
00359 Psi.MatvecLT(h, g, res);
00360 for(i = 0; i <= N; ++i) {
00361 res[i] -= 0.5 * h *
00362 (Psi.GetM().elem(i, 0) * g[0] + Psi.GetM().elem(i, i) * g[i]);
00363 }
00364 }
00365
00366 void SmoluchowskiCalcFast::calc_L4(const double *g, double *res) const
00367 {
00368 unsigned i;
00369 Psi.MatvecUT(h, g, res, Skeleton::transpose);
00370 for(i = 0; i <= N; ++i) {
00371 res[i] -= 0.5 * h *
00372 (Psi.GetM().elem(i, i) * g[i] + Psi.GetM().elem(N, i) * g[N]);
00373 }
00374 }
00375
00376 void SmoluchowskiCalcFast::calc_L5(const double *g, double *res) const
00377 {
00378 unsigned i;
00379 K.MatvecF(h, g, res, tmp_mv);
00380 for(i = 0; i <= N; ++i) {
00381 res[i] -= 0.5 * h *
00382 (K.GetM().elem(i, 0) * g[0] + K.GetM().elem(i, N) * g[N]);
00383 }
00384 }
00385
00386
00387
00388
00389 SmoluchowskiOperator::SmoluchowskiOperator(SmoluchowskiCalc& calc)
00390 : smol_base(calc)
00391 {
00392 N = smol_base.GetN();
00393 H = smol_base.GetH();
00394 h = H / double(N);
00395 tmp_l = new double[N + 1];
00396 }
00397
00398 SmoluchowskiOperator::~SmoluchowskiOperator()
00399 {
00400 delete [] tmp_l;
00401 }
00402
00403
00404
00405 SmoluchowskiNonLinearOperator::SmoluchowskiNonLinearOperator(
00406 SmoluchowskiCalc& calc)
00407 : SmoluchowskiOperator(calc)
00408 { }
00409
00410 SmoluchowskiNonLinearOperator::~SmoluchowskiNonLinearOperator()
00411 { }
00412
00413 void SmoluchowskiNonLinearOperator::Apply(const double *g, double *res)
00414 {
00415 unsigned i;
00416 for(i = 0; i <= N; ++i) {
00417 res[i] = double(i) * h;
00418 }
00419 smol_base.calc_L3(res, tmp_l);
00420 smol_base.calc_L4(g, res);
00421 for(i = 1; i <= N; ++i) {
00422 res[i] -= g[i] / (double(i) * h) * tmp_l[i];
00423 }
00424 smol_base.calc_L5(g, tmp_l);
00425 for(i = 0; i <= N; ++i) {
00426 res[i] -= g[i] * tmp_l[i];
00427 }
00428 smol_base.calc_L1(g, g, tmp_l);
00429 for(i = 0; i <= N; ++i) {
00430 res[i] = -(res[i] + 0.5 * tmp_l[i]);
00431 }
00432 }
00433
00434
00435
00436 SmoluchowskiLinearOperator::SmoluchowskiLinearOperator(
00437 SmoluchowskiCalc& calc, const double *st)
00438 : SmoluchowskiOperator(calc)
00439 {
00440 unsigned i;
00441 c0 = new double[N + 1];
00442 for(i = 0; i <= N; ++i) {
00443 c0[i] = st[i];
00444 }
00445 smol_base.fix_f_arg(c0);
00446 a = new double[N + 1];
00447 calc_a();
00448 }
00449
00450 SmoluchowskiLinearOperator::~SmoluchowskiLinearOperator()
00451 {
00452 delete [] c0;
00453 delete [] a;
00454 }
00455
00456 const double *SmoluchowskiLinearOperator::Get_c0() const
00457 {
00458 return c0;
00459 }
00460
00461 const double *SmoluchowskiLinearOperator::Get_a() const
00462 {
00463 return a;
00464 }
00465
00466
00467 void SmoluchowskiLinearOperator::calc_a()
00468 {
00469 unsigned i;
00470 for(i = 0; i <= N; ++i) {
00471 a[i] = double(i) * h;
00472 }
00473 smol_base.calc_L3(a, tmp_l);
00474 smol_base.calc_L5(c0, a);
00475 for(i = 1; i <= N; ++i) {
00476 a[i] += tmp_l[i] / ((double)(i) * h);
00477 }
00478 }
00479
00480 void SmoluchowskiLinearOperator::Apply(const double *g, double *res)
00481 {
00482 unsigned i;
00483 smol_base.calc_L5(g, res);
00484 smol_base.calc_L1_fix_f(g, tmp_l);
00485 for(i = 0; i <= N; ++i) {
00486 res[i] = c0[i] * res[i] - tmp_l[i];
00487 }
00488 smol_base.calc_L4(g, tmp_l);
00489 for(i = 0; i <= N; ++i) {
00490 res[i] += a[i] * g[i] - tmp_l[i];
00491 }
00492 }
00493
00494 void SmoluchowskiLinearOperator::ApplyAdjoint(const double *g, double *res)
00495 {
00496 unsigned i;
00497 for(i = 0; i <= N; ++i) {
00498 tmp_l[i] = c0[i] * g[i];
00499 }
00500 smol_base.calc_L5(tmp_l, res);
00501 smol_base.calc_L2_fix_f(g, tmp_l);
00502 for(i = 0; i <= N; ++i) {
00503 res[i] -= tmp_l[i];
00504 }
00505 smol_base.calc_L3(g, tmp_l);
00506 for(i = 0; i <= N; ++i) {
00507 res[i] += a[i] * g[i] - tmp_l[i];
00508 }
00509 }