00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040 #include "matrixbk.hh"
00041 #include <strstream>
00042 #include <algorithm>
00043
00044 #ifdef __DEBUGMAT1
00045 #include <fstream>
00046 #endif
00047
00048
00049 #include <cstdlib>
00050
00051
00052
00053
00054
00055 extern bk_messages matERROR;
00056 extern bk_messages matDEBUG;
00057
00058
00059
00060
00061
00062
00063
00064 void four1(
00065 complr4 data[],
00066 int32 fftlength,
00067 int32 isign);
00068
00069
00070
00071 #ifdef __USE_VECLIB_LIBRARY__
00072 extern "C" { int32 sgemm (char* , char*, int32*, int32*, int32*,
00073 real4*, real4*, int32*,
00074 real4*, int32*, real4*,
00075 real4*, int32*, int32, int32); }
00076 extern "C" { int32 dgemm (char* , char*, int32*, int32*, int32*,
00077 real8*, real8*, int32*,
00078 real8*, int32*, real8*,
00079 real8*, int32*, int32, int32); }
00080 extern "C" { int32 cgemm (char* , char*, int32*, int32*, int32*,
00081 complr4*, complr4*, int32*,
00082 complr4*, int32*, complr4*,
00083 complr4*, int32*, int32, int32); }
00084 extern "C" { int32 zgemm (char* , char*, int32*, int32*, int32*,
00085 complr8*, complr8*, int32*,
00086 complr8*, int32*, complr8*,
00087 complr8*, int32*, int32, int32); }
00088 #endif
00089
00090 #ifdef __USE_VECLIB_LIBRARY__
00091 extern "C" { int32 c1dfft (complr4*, int32*, real4*, int32*, int32*); }
00092 extern "C" { int32 c2dfft (complr4*, int32*, int32*, int32*, int32*, int32*); }
00093
00094 #endif
00095
00096
00097
00098 #ifdef __USE_FFTW_LIBRARY__
00099 #include <fftw3.h>
00100 #endif
00101
00102 #ifdef __USE_LAPACK_LIBRARY__
00103 extern "C" { int32 spotrf (char*, int32*, real4*, int32*, int32*, int32); }
00104 extern "C" { int32 spotri (char*, int32*, real4*, int32*, int32*, int32); }
00105 extern "C" { int32 spotrs (char*, int32*, int32*, real4*, int32*,
00106 real4*, int32*, int32*, int32); }
00107 extern "C" { int32 dpotrf (char*, int32*, real8*, int32*, int32*, int32); }
00108 extern "C" { int32 dpotri (char*, int32*, real8*, int32*, int32*, int32); }
00109 extern "C" { int32 dpotrs (char*, int32*, int32*, real8*, int32*,
00110 real8*, int32*, int32*, int32); }
00111 #endif
00112
00113
00114
00115
00116
00117
00118
00119
00120 inline bool myispower2(uint w)
00121 {if (w==2 || w==4 || w==8 || w==16 ||
00122 w==32 || w==64 || w==128 || w==256 ||
00123 w==512 || w==1024 || w==2048 || w==4096 ||
00124 w==8192 || w==16384 || w==32768 )
00125 return true;
00126 return false;}
00127
00128
00129
00130
00131
00132
00133
00134
00135 void matassert(
00136 const ifstream &str,
00137 const char* ifilename,
00138 const char* callingfile,
00139 int32 linenumber)
00140 {
00141 if (!str)
00142 {
00143 matERROR << "Cannot open file: \"" << ifilename
00144 << "\" (input). Source: \"" << callingfile
00145 << "\", line: " << linenumber;
00146 matERROR.print();
00147 }
00148 #ifdef __DEBUGMAT2
00149 matDEBUG << "OK input stream associated with file: \""
00150 << ifilename << "\"";
00151 matDEBUG.print();
00152 #endif
00153 }
00154
00155
00156
00157
00158
00159 void matassert(
00160 const ofstream &str,
00161 const char* ofilename,
00162 const char* callingfile,
00163 int32 linenumber)
00164 {
00165 if (!str)
00166 {
00167 matERROR << "Cannot open file: \"" << ofilename
00168 << "\" (output). Source: \"" << callingfile
00169 << "\", line: " << linenumber
00170 << "; OVERWRIT OFF/non existing directory?";
00171 matERROR.print();
00172 }
00173 #ifdef __DEBUGMAT2
00174 matDEBUG << "OK output stream associated with file: \""
00175 << ofilename << "\"";
00176 matDEBUG.print();
00177 #endif
00178 }
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192 void fileci2tomatcr4(
00193 matrix<complr4> &Result,
00194 const char *file,
00195 uint filelines,
00196 window win,
00197 window winoffset)
00198 {
00199 #ifdef __DEBUGMAT2
00200 matDEBUG.print("ROUTINE: fileci2tomatcr4.");
00201 #endif
00202
00203
00204 win.linelo = win.linelo - winoffset.linelo + 1;
00205 win.linehi = win.linehi - winoffset.linelo + 1;
00206 win.pixlo = win.pixlo - winoffset.pixlo + 1;
00207 win.pixhi = win.pixhi - winoffset.pixlo + 1;
00208
00209
00210
00211
00212
00213 #ifdef __NO_IOS_BINARY__
00214 ifstream ifile(file, ios::in);
00215 #else
00216 ifstream ifile(file, ios::in | ios::binary);
00217 #endif
00218
00219
00220 ifile.seekg(0,ios::end);
00221 matassert(ifile,file,__FILE__,__LINE__);
00222 const uint sizefile = ifile.tellg();
00223
00224 if (sizefile==0)
00225 matERROR.print("filesize==0...");
00226
00227
00228
00229 const uint pixelsxbytes = sizefile/filelines;
00230 const uint SIZE = sizeof(compli16);
00231 const uint filepixels = pixelsxbytes/SIZE;
00232 const uint sizepixel = pixelsxbytes/filepixels;
00233 const uint lines = win.lines();
00234 const uint pixels = win.pixels();
00235
00236
00237 const uint start = ((win.linelo-1)*filepixels+win.pixlo-1)*sizepixel;
00238
00239 #ifdef __DEBUGMAT2
00240 if (win.linelo<=0)
00241 matERROR.print("minimum line <= 0, (should be 1?).");
00242 if (win.pixlo<=0)
00243 matERROR.print("minimum pixel <= 0, (should be 1?).");
00244 if (win.linelo>win.linehi)
00245 matERROR.print("minimum line > max line.");
00246 if (win.pixlo>win.pixhi)
00247 matERROR.print("minimum pixel > max pixel.");
00248 if (win.linehi>filelines)
00249 matERROR.print("max. line is larger than on file.");
00250 if (win.pixhi>filepixels)
00251 matERROR.print("max. pixel is larger than on file.");
00252 if (SIZE != sizepixel)
00253 matERROR.print("Formatflag not consistent with file.");
00254 #endif
00255
00256 #ifdef __DEBUGMAT2
00257 if (Result.lines() < lines || Result.pixels() < pixels)
00258 matERROR.print("code ?::fileci2tomatcr4: Result matrix not ok.");
00259 if (Result.lines() != lines || Result.pixels() != pixels)
00260 matDEBUG.print("debug info: fileci2tomatcr4: matrix not fully filled.");
00261 #endif
00262
00263 int16 el[2];
00264
00265 for (int32 lin=0; lin<lines; ++lin)
00266 {
00267 ifile.seekg(start+filepixels*lin*sizepixel,ios::beg);
00268 for (int32 pix=0; pix<pixels; ++pix)
00269 {
00270 ifile.read((char*)el,SIZE);
00271 Result(lin,pix)=complr4(el[0],el[1]);
00272
00273 }
00274 }
00275 ifile.close();
00276 }
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286 matrix<real4> magnitude(
00287 const matrix<complr4> &A)
00288 {
00289 #ifdef __DEBUGMAT2
00290 matDEBUG.print("magnitude.");
00291 #endif
00292
00293 matrix<real4> Result(A.lines(),A.pixels());
00294 complr4 *pntA = A[0];
00295 real4 *pntR = Result[0];
00296 for (register int32 i=0; i<Result.size(); i++)
00297 (*pntR++) = sqrt(sqr((*pntA).real()) + sqr((*pntA++).imag()));
00298 return Result;
00299 }
00300
00301
00302
00303
00304
00305
00306
00307
00308 matrix<real4> intensity(
00309 const matrix<complr4> &A)
00310 {
00311 #ifdef __DEBUGMAT2
00312 matDEBUG.print("intensity.");
00313 #endif
00314
00315 matrix<real4> Result(A.lines(),A.pixels());
00316 complr4 *pntA = A[0];
00317 real4 *pntR = Result[0];
00318 for (register int32 i=0; i<Result.size(); i++)
00319 (*pntR++) = sqr((*pntA).real()) + sqr((*pntA++).imag());
00320 return Result;
00321 }
00322
00323
00324
00325 #ifdef __USE_LAPACK_LIBRARY__
00326
00327
00328
00329
00330
00331
00332 void choles(
00333 matrix<real4> &A)
00334 {
00335 #ifdef __DEBUGMAT2
00336 matDEBUG.print("choles: LAPACK");
00337 #endif
00338 #ifdef __DEBUGMAT1
00339 if (A.lines() != A.pixels())
00340 matERROR.print("only symmetrical matrixes for cholesky.");
00341 #endif
00342
00343 int32 lda= A.pixels();
00344 int32 n = lda;
00345 int32 ierr;
00346 spotrf("U",&n,A[0],&lda,&ierr,1);
00347
00348 #ifdef __DEBUGMAT1
00349 if (ierr == 0)
00350 matDEBUG.print("choles: ok. lower contains factorisation.");
00351 else if (ierr < 0)
00352 matERROR.print("choles: the kth element is not ok.");
00353 else
00354 matERROR.print("choles: not positive definite.");
00355 #endif
00356 }
00357
00358
00359
00360
00361 #else
00362
00363
00364
00365
00366
00367
00368
00369
00370 void choles(matrix<real4> &A)
00371 {
00372 #ifdef __DEBUGMAT2
00373 matDEBUG.print("choles: internal");
00374 #endif
00375 const int32 N = A.lines();
00376 register real4 sum;
00377 for (register int32 i=0; i<N; ++i)
00378 {
00379 for (register int32 j=i; j<N; ++j)
00380 {
00381 sum = A[i][j];
00382 for (register int32 k=i-1; k>=0; --k)
00383 {
00384 sum -= A[i][k] * A[j][k];
00385 }
00386 if (i == j)
00387 {
00388 if (sum <= 0.) {matERROR.print("choles: internal: A not pos. def.");}
00389 A[i][i] = sqrt(sum);
00390 }
00391 else
00392 {
00393 A[j][i] = sum/A[i][i];
00394 }
00395 }
00396 }
00397 }
00398 #endif
00399
00400
00401 #ifdef __USE_LAPACK_LIBRARY__
00402
00403
00404
00405
00406
00407
00408 void choles(
00409 matrix<real8> &A)
00410 {
00411 #ifdef __DEBUGMAT2
00412 matDEBUG.print("choles: LAPACK");
00413 #endif
00414 #ifdef __DEBUGMAT1
00415 if (A.lines() != A.pixels())
00416 matERROR.print("only symmetrical matrixes for cholesky.");
00417 #endif
00418
00419 int32 lda= A.pixels();
00420 int32 n = lda;
00421 int32 ierr;
00422 dpotrf("U",&n,A[0],&lda,&ierr,1);
00423
00424 #ifdef __DEBUGMAT1
00425 if (ierr == 0)
00426 matDEBUG.print("choles: ok. lower contains factorisation.");
00427 else if (ierr < 0)
00428 matERROR.print("choles: the kth element is not ok.");
00429 else
00430 matERROR.print("choles: not positive definite.");
00431 #endif
00432 }
00433
00434
00435
00436 #else
00437
00438
00439
00440
00441
00442
00443
00444
00445 void choles(matrix<real8> &A)
00446 {
00447 #ifdef __DEBUGMAT2
00448 matDEBUG.print("choles: internal");
00449 #endif
00450 const int32 N = A.lines();
00451 register real8 sum;
00452 for (register int32 i=0; i<N; ++i)
00453 {
00454 for (register int32 j=i; j<N; ++j)
00455 {
00456 sum = A[i][j];
00457 for (register int32 k=i-1; k>=0; --k)
00458 {
00459 sum -= A[i][k] * A[j][k];
00460 }
00461 if (i == j)
00462 {
00463 if (sum <= 0.) {matERROR.print("choles: internal: A not pos. def.");}
00464 A[i][i] = sqrt(sum);
00465 }
00466 else
00467 {
00468 A[j][i] = sum/A[i][i];
00469 }
00470 }
00471 }
00472 }
00473 #endif
00474
00475
00476 #ifdef __USE_LAPACK_LIBRARY__
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487 void solvechol(
00488 const matrix<real4> &A,
00489 matrix<real4> &B)
00490 {
00491 #ifdef __DEBUGMAT2
00492 matDEBUG.print("solvechol::LAPACK");
00493 #endif
00494 #ifdef __DEBUGMAT1
00495 if (A.lines() != A.pixels())
00496 matERROR.print("only symmetrical matrices for cholesky.");
00497 if (A.lines() != B.lines())
00498 matERROR.print("solvechol: must same size a,b.");
00499 if (B.pixels() != 1)
00500 matERROR.print("solvechol: NOT possible simultaneous solution in this way.");
00501 #endif
00502
00503 int32 lda= A.pixels();
00504 int32 ldb= B.lines();
00505 int32 nrhs= B.pixels();
00506 int32 n = lda;
00507 int32 ierr;
00508 spotrs("U",&n,&nrhs,A[0],&lda,B[0],&ldb,&ierr,1);
00509
00510 #ifdef __DEBUGMAT1
00511 if (ierr == 0)
00512 matDEBUG.print("solvechol: ok, rhs contains solution.");
00513 else if (ierr < 0)
00514 matERROR.print("solvechol: the kth element is not ok.");
00515 else
00516 matERROR.print("solvechol: impossible ierr.");
00517 #endif
00518 }
00519
00520
00521
00522 #else
00523
00524
00525
00526
00527
00528
00529
00530
00531 void solvechol(
00532 const matrix<real4> &A,
00533 matrix<real4> &B)
00534 {
00535 #ifdef __DEBUGMAT2
00536 matDEBUG.print("solvechol::INTERNAL");
00537 #endif
00538 #ifdef __DEBUGMAT1
00539 if (A.lines() != A.pixels())
00540 matERROR.print("solvechol: INTERNAL: symmetrical matrices for cholesky.");
00541 if (A.lines() != B.lines())
00542 matERROR.print("solvechol: must same size a,b.");
00543 if (B.pixels() != 1)
00544 matERROR.print("solvechol: NOT possible simultaneous solution in this way.");
00545 #endif
00546 const int32 N = A.lines();
00547 register real4 sum;
00548 register int32 i,j;
00549
00550 for (i=0; i<N; ++i)
00551 {
00552 sum = B[i][0];
00553 for (j=i-1; j>=0; --j)
00554 {
00555 sum -= A[i][j]*B[j][0];
00556 }
00557 B[i][0] = sum/A[i][i];
00558 }
00559
00560 for (i=N-1; i>=0; --i)
00561 {
00562 sum = B[i][0];
00563 for (j=i+1; j<N; ++j)
00564 {
00565 sum -= A[j][i]*B[j][0];
00566 }
00567 B[i][0] = sum/A[i][i];
00568 }
00569 }
00570 #endif
00571
00572
00573 #ifdef __USE_LAPACK_LIBRARY__
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584 void solvechol(
00585 const matrix<real8> &A,
00586 matrix<real8> &B)
00587 {
00588 #ifdef __DEBUGMAT2
00589 matDEBUG.print("solvechol: LAPACK: there was problem here, ldb,nrhs.");
00590 #endif
00591 #ifdef __DEBUGMAT1
00592 if (A.lines() != A.pixels())
00593 matERROR.print("only symmetrical matrixes for cholesky.");
00594 if (A.lines() != B.lines())
00595 matERROR.print("solvechol: must same size A,B.");
00596 if (B.pixels() != 1)
00597 matERROR.print("solvechol: NOT possible simultaneous solution in this way.");
00598 #endif
00599
00600 int32 lda = A.pixels();
00601 int32 ldb = B.lines();
00602 int32 nrhs = B.pixels();
00603 int32 n = lda;
00604 int32 ierr;
00605 dpotrs("U",&n,&nrhs,A[0],&lda,B[0],&ldb,&ierr,1);
00606
00607 #ifdef __DEBUGMAT1
00608 if (ierr == 0)
00609 matDEBUG.print("solvechol: ok, rhs contains solution.");
00610 else if (ierr < 0)
00611 matERROR.print("solvechol: the kth element is not ok.");
00612 else
00613 matERROR.print("solvechol: impossible ierr.");
00614 #endif
00615 }
00616
00617
00618
00619 #else
00620
00621
00622
00623
00624
00625
00626
00627
00628 void solvechol(
00629 const matrix<real8> &A,
00630 matrix<real8> &B)
00631 {
00632 #ifdef __DEBUGMAT2
00633 matDEBUG.print("solvechol::INTERNAL");
00634 #endif
00635 #ifdef __DEBUGMAT1
00636 if (A.lines() != A.pixels())
00637 matERROR.print("solvechol: INTERNAL: symmetrical matrices for cholesky.");
00638 if (A.lines() != B.lines())
00639 matERROR.print("solvechol: must same size a,b.");
00640 if (B.pixels() != 1)
00641 matERROR.print("solvechol: NOT possible simultaneous solution in this way.");
00642 #endif
00643 const int32 N = A.lines();
00644 register real8 sum;
00645 register int32 i,j;
00646
00647 for (i=0; i<N; ++i)
00648 {
00649 sum = B[i][0];
00650 for (j=i-1; j>=0; --j)
00651 {
00652 sum -= A[i][j]*B[j][0];
00653 }
00654 B[i][0] = sum/A[i][i];
00655 }
00656
00657 for (i=N-1; i>=0; --i)
00658 {
00659 sum = B[i][0];
00660 for (j=i+1; j<N; ++j)
00661 {
00662 sum -= A[j][i]*B[j][0];
00663 }
00664 B[i][0] = sum/A[i][i];
00665 }
00666 }
00667 #endif
00668
00669
00670 #ifdef __USE_LAPACK_LIBRARY__
00671
00672
00673
00674
00675
00676
00677 void invertchol(
00678 matrix<real4> &A)
00679 {
00680 #ifdef __DEBUGMAT2
00681 matDEBUG.print("invertchol: LAPACK");
00682 #endif
00683 #ifdef __DEBUGMAT1
00684 if (A.lines() != A.pixels())
00685 matERROR.print("only symmetrical matrixes for cholesky.");
00686 #endif
00687
00688 int32 lda= A.pixels();
00689 int32 n = lda;
00690 int32 ierr;
00691 spotri("U",&n,A[0],&lda,&ierr,1);
00692
00693 #ifdef __DEBUGMAT1
00694 if (ierr == 0)
00695 matDEBUG.print("invertcholes: ok, lower contains inverse.");
00696 else if (ierr < 0)
00697 matERROR.print("invertcholes: the kth element is not ok.");
00698 else
00699 matERROR.print("invertcholes: singular.");
00700 #endif
00701 }
00702
00703
00704
00705
00706 #else
00707
00708
00709
00710
00711
00712
00713
00714 void invertchol(
00715 matrix<real4> &A)
00716 {
00717 const int32 N = A.lines();
00718 register real4 sum;
00719 register int32 i,j,k;
00720
00721 for (i=0; i<N; ++i)
00722 {
00723 A[i][i] = 1./A[i][i];
00724 for (j=i+1; j<N; ++j)
00725 {
00726 sum = 0.;
00727 for (k=i; k<j; ++k)
00728 {
00729 sum -= A[j][k]*A[k][i];
00730 }
00731 A[j][i]=sum/A[j][j];
00732 }
00733 }
00734
00735 for (i=0; i<N; ++i)
00736 {
00737 for (j=i; j<N; ++j)
00738 {
00739 sum = 0.;
00740 for (k=j; k<N; ++k)
00741 {
00742 sum += A[k][i]*A[k][j];
00743 }
00744 A[j][i] = sum;
00745 }
00746 }
00747 }
00748 #endif
00749
00750
00751 #ifdef __USE_LAPACK_LIBRARY__
00752
00753
00754
00755
00756
00757
00758 void invertchol(
00759 matrix<real8> &A)
00760 {
00761 #ifdef __DEBUGMAT2
00762 matDEBUG.print("invertchol: LAPACK");
00763 #endif
00764 #ifdef __DEBUGMAT1
00765 if (A.lines() != A.pixels())
00766 matERROR.print("only symmetrical matrixes for cholesky.");
00767 #endif
00768
00769 int32 lda= A.pixels();
00770 int32 n = lda;
00771 int32 ierr;
00772 dpotri("U",&n,A[0],&lda,&ierr,1);
00773
00774 #ifdef __DEBUGMAT1
00775 if (ierr == 0)
00776 matDEBUG.print("invertcholes: ok, lower contains inverse.");
00777 else if (ierr < 0)
00778 matERROR.print("invertcholes: the kth element is not ok.");
00779 else
00780 matERROR.print("invertcholes: singular.");
00781 #endif
00782 }
00783
00784
00785
00786 #else
00787
00788
00789
00790
00791
00792
00793
00794 void invertchol(
00795 matrix<real8> &A)
00796 {
00797 const int32 N = A.lines();
00798 register real8 sum;
00799 register int32 i,j,k;
00800
00801 for (i=0; i<N; ++i)
00802 {
00803 A[i][i] = 1./A[i][i];
00804 for (j=i+1; j<N; ++j)
00805 {
00806 sum = 0.;
00807 for (k=i; k<j; ++k)
00808 {
00809 sum -= A[j][k]*A[k][i];
00810 }
00811 A[j][i]=sum/A[j][j];
00812 }
00813 }
00814
00815 for (i=0; i<N; ++i)
00816 {
00817 for (j=i; j<N; ++j)
00818 {
00819 sum = 0.;
00820 for (k=j; k<N; ++k)
00821 {
00822 sum += A[k][i]*A[k][j];
00823 }
00824 A[j][i] = sum;
00825 }
00826 }
00827 }
00828 #endif
00829
00830
00831
00832
00833
00834
00835
00836 matrix<complr4> norm(
00837 const matrix<complr4> &A)
00838 {
00839 #ifdef __DEBUGMAT2
00840 matDEBUG.print("norm");
00841 #endif
00842
00843 matrix<complr4> Result(A.lines(),A.pixels());
00844 complr4 *pntA = A[0];
00845 complr4 *pntR = Result[0];
00846 for (register int32 i=0; i<Result.size(); i++)
00847 (*pntR++) = complr4(norm(*pntA++));
00848 return Result;
00849 }
00850
00851
00852
00853
00854
00855
00856
00857 real4 norm2(
00858 const matrix<complr4> &A)
00859 {
00860 #ifdef __DEBUGMAT2
00861 matDEBUG.print("norm2");
00862 #endif
00863
00864 real4 n=0.;
00865 complr4 *pntA = A[0];
00866 for (register int32 i=0; i<A.size(); i++)
00867 n += norm(*pntA++);
00868 return n;
00869 }
00870
00871
00872
00873
00874
00875
00876
00877 real4 norm2(
00878 const matrix<real4> &A)
00879 {
00880 #ifdef __DEBUGMAT2
00881 matDEBUG.print("norm2");
00882 #endif
00883
00884 real4 n=0.;
00885 real4 *pntA = A[0];
00886 for (register int32 i=0; i<A.size(); i++)
00887 n += sqr(*pntA++);
00888 return n;
00889 }
00890
00891
00892
00893
00894
00895
00896
00897 matrix<real4> abs(
00898 const matrix<real4> &A)
00899 {
00900 #ifdef __DEBUGMAT2
00901 matDEBUG.print("abs");
00902 #endif
00903 matrix<real4> Result(A.lines(),A.pixels());
00904
00905 real4 *pntA = A[0];
00906 real4 *pntR = Result[0];
00907 for (register int32 i=0; i<Result.size(); i++)
00908 (*pntR++) = (abs((*pntA++)));
00909 return Result;
00910 }
00911
00912
00913
00914
00915
00916
00917
00918 matrix<real8> abs(
00919 const matrix<real8> &A)
00920 {
00921 #ifdef __DEBUGMAT2
00922 matDEBUG.print("abs");
00923 #endif
00924 matrix<real8> Result(A.lines(),A.pixels());
00925
00926 real8 *pntA = A[0];
00927 real8 *pntR = Result[0];
00928 for (register int32 i=0; i<Result.size(); i++)
00929 (*pntR++) = (abs((*pntA++)));
00930 return Result;
00931 }
00932
00933
00934
00935
00936
00937
00938
00939 matrix<real4> real(
00940 const matrix<complr4> &A)
00941 {
00942 matrix<real4> Result(A.lines(),A.pixels());
00943 real4 *pntR = Result[0];
00944 complr4 *pntA = A[0];
00945 #ifdef __DEBUGMAT2
00946 matDEBUG.print("real. do not initialize memory (no alloc&and init).");
00947 #endif
00948
00949 for (register int32 i=0; i<A.size(); i++)
00950 (*pntR++) = real(*pntA++);
00951 return Result;
00952 }
00953
00954
00955
00956
00957
00958
00959
00960 matrix<real4> imag(
00961 const matrix<complr4> &A)
00962 {
00963 matrix<real4> Result(A.lines(),A.pixels());
00964 real4 *pntR = Result[0];
00965 complr4 *pntA = A[0];
00966 #ifdef __DEBUGMAT2
00967 matDEBUG.print("imag. do not initialize memory (no alloc).");
00968 #endif
00969
00970 for (register int32 i=0; i<A.size(); i++)
00971 (*pntR++) = imag(*pntA++);
00972 return Result;
00973 }
00974
00975
00976
00977
00978
00979
00980
00981 matrix<real4> angle(
00982 const matrix<complr4> &A)
00983 {
00984 matrix<real4> Result(A.lines(),A.pixels());
00985 real4 *pntR = Result[0];
00986 complr4 *pntA = A[0];
00987 #ifdef __DEBUGMAT2
00988 matDEBUG.print("angle. do not initialize memory (no alloc).");
00989 #endif
00990
00991
00992 for (register int32 i=0; i<A.size(); i++)
00993 (*pntR++) = arg(*pntA++);
00994
00995 return Result;
00996 }
00997
00998
00999
01000
01001
01002
01003
01004 matrix<complr4> angle2cmplx(
01005 const matrix<real4> &A)
01006 {
01007 matrix<complr4> Result(A.lines(),A.pixels());
01008 complr4 *pntR = Result[0];
01009 real4 *pntA = A[0];
01010 #ifdef __DEBUGMAT2
01011 matDEBUG.print("angle2cmplx. do not initialize memory (no alloc).");
01012 #endif
01013
01014 for (register int32 i=0; i<A.size(); i++)
01015 (*pntR++) = complr4(cos(*pntA),sin(*pntA++));
01016 return Result;
01017 }
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027 void dotmultconjphase(
01028 matrix<complr4> &cint,
01029 const matrix<real4> &refpha)
01030 {
01031 complr4 *pntC = cint[0];
01032 real4 *pntA = refpha[0];
01033 #ifdef __DEBUGMAT2
01034 matDEBUG.print("dotmultconjphase.");
01035 #endif
01036 #ifdef __DEBUGMAT1
01037 if (cint.lines() != refpha.lines())
01038 matERROR.print("dotmultconjphase: not same number of lines.");
01039 if (cint.pixels() != refpha.pixels())
01040 matERROR.print("dotmultconjphase: not same number of pixels.");
01041 #endif
01042
01043
01044
01045
01046
01047
01048 for (register int32 i=0; i<cint.size(); i++)
01049 (*pntC++) *= complr4(cos(*pntA),-sin(*pntA++));
01050 }
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061 matrix<complr4> coherence(
01062 const matrix<complr4> &CINT,
01063 const matrix<complr4> &NORMS,
01064 uint winL,
01065 uint winP)
01066 {
01067
01068
01069 #define PLACE2 // fastest
01070
01071
01072 #ifdef __DEBUGMAT2
01073 matDEBUG.print("coherence");
01074 if (!winL>=winP)
01075 matDEBUG.print("coherence: estimator window size L<P not very efficiently programmed.");
01076 #endif
01077 #ifdef __DEBUGMAT1
01078 if (CINT.lines() != NORMS.lines() ||
01079 CINT.pixels() != NORMS.pixels() )
01080 matERROR.print("not same dimensions.");
01081 #endif
01082
01083
01084
01085 matrix<complr4> Result(CINT.lines()-winL+1,CINT.pixels());
01086 int32 leadingzeros = (winP-1)/2;
01087 int32 trailingzeros = (winP)/2;
01088
01089
01090 #ifdef PLACE
01091 complr4 sum = 0.;
01092 complr4 power = 0.;
01093 for (register int32 i=0; i<Result.lines(); i++)
01094 {
01095 for (register int32 j=leadingzeros; j<Result.pixels()-trailingzeros; j++)
01096 {
01097 for (register int32 k=i; k<i+winL; k++)
01098 {
01099 for (register int32 l=j-leadingzeros; l<j-leadingzeros+winP; l++)
01100 {
01101 sum += CINT(k,l);
01102 power += NORMS(k,l);
01103 }
01104 }
01105 real4 p = power.real() * power.imag();
01106 Result(i,j) = (p > 0.0) ? sum/sqrt(p) : 0.0;
01107 sum = 0.;
01108 power = 0.;
01109 }
01110 }
01111 #undef PLACE
01112 #endif
01113
01114
01115
01116
01117
01118 #ifdef PLACE2
01119 register int32 i,j,k,l;
01120 register complr4 sum;
01121 register complr4 power;
01122 for (j=leadingzeros; j<Result.pixels()-trailingzeros; j++)
01123 {
01124 sum = complr4(0.);
01125 power = complr4(0.);
01126
01127 for (k=0; k<winL; k++)
01128 {
01129 for (l=j-leadingzeros; l<j-leadingzeros+winP; l++)
01130 {
01131 sum += CINT(k,l);
01132 power += NORMS(k,l);
01133 }
01134 }
01135 real4 p = power.real() * power.imag();
01136 Result(0,j) = (p > 0.0) ? sum/sqrt(p) : 0.0;
01137
01138
01139 for (i=0; i<Result.lines()-1; i++)
01140 {
01141 for (l=j-leadingzeros; l<j-leadingzeros+winP; l++)
01142 {
01143 sum += (CINT(i+winL,l) - CINT(i,l));
01144 power += (NORMS(i+winL,l) - NORMS(i,l));
01145 }
01146 real4 p = power.real() * power.imag();
01147 Result(i+1,j) = (p > 0.0) ? sum/sqrt(p) : 0.0;
01148 }
01149 }
01150 #undef PLACE2
01151 #endif
01152
01153
01154
01155
01156
01157
01158 #ifdef FFTS
01159 matDEBUG.print("THIS SEEMS TO GO WRONG SOMEWHERE, writes all zeros, computes ok?");
01160 bool newline = false;
01161 bool nonewline = false;
01162
01163 int32 SIZEBLOCKL = 256;
01164 int32 SIZEBLOCKP = 1024;
01165 int32 sizeL = SIZEBLOCKL;
01166 int32 sizeP = SIZEBLOCKP;
01167 if (CINT.lines() < SIZEBLOCKL)
01168 {
01169 SIZEBLOCKL = nextpow2(CINT.lines());
01170 sizeL = CINT.lines();
01171 nonewline = true;
01172 }
01173 if (CINT.pixels() < SIZEBLOCKP)
01174 {
01175 SIZEBLOCKP = nextpow2(CINT.pixels());
01176 sizeP = CINT.pixels();
01177 }
01178
01179 const int32 twoL = 2*SIZEBLOCKL;
01180 const int32 twoP = 2*SIZEBLOCKP;
01181 const complr4 cr4one = complr4(1.,0);
01182
01183 matrix<complr4> CINT2(twoL, twoP);
01184 matrix<complr4> NORMS2(twoL, twoP);
01185 matrix<complr4> BLOCK(twoL, twoP);
01186 register int32 i,j;
01187 for (i=0; i<winL; i++)
01188 for (j=0; j<winP; j++)
01189 BLOCK(i,j) = cr4one;
01190 fft2d(BLOCK);
01191 BLOCK = conj(BLOCK);
01192
01193 window win = {0, sizeL-1, 0, sizeP-1};
01194 window win2 = {0, sizeL-1, 0, sizeP-1};
01195 register int32 i2 = 0;
01196 register int32 j2 = leadingzeros;
01197
01198 for (;;)
01199 {
01200 #ifdef __DEBUGMAT2
01201 cout << "time for coherence init:\n";
01202 printcpu();
01203 #endif
01204
01205 cout << "BERT: win:"
01206 << win.linelo << " "
01207 << win.linehi << " "
01208 << win.pixlo << " "
01209 << win.pixhi << "\n";
01210 cout << "BERT: win2:"
01211 << win2.linelo << " "
01212 << win2.linehi << " "
01213 << win2.pixlo << " "
01214 << win2.pixhi << "\n";
01215 cout << "BERT: i2,j2 sizeL,p:"
01216 << i2 << " "
01217 << j2 << " "
01218 << sizeL << " "
01219 << sizeP << "\n";
01220
01221
01222 CINT2.setdata(win2,CINT,win);
01223 NORMS2.setdata(win2,NORMS,win);
01224 #ifdef __DEBUGMAT2
01225 cout << "time for setdata:\n";
01226 printcpu();
01227 #endif
01228 fft2d(CINT2);
01229 fft2d(NORMS2);
01230 #ifdef __DEBUGMAT2
01231 cout << "time for fft:\n";
01232 printcpu();
01233 #endif
01234
01235
01236 for (i=0; i<CINT2.lines(); i++)
01237 {
01238 for (j=0; j<CINT2.pixels(); j++)
01239 {
01240 CINT2(i,j) *= BLOCK(i,j);
01241 NORMS2(i,j) *= BLOCK(i,j);
01242 }
01243 }
01244
01245 #ifdef __DEBUGMAT2
01246 cout << "time for convol.:\n";
01247 printcpu();
01248 #endif
01249 ifft2d(CINT2);
01250 ifft2d(NORMS2);
01251 #ifdef __DEBUGMAT2
01252 cout << "time for ifft.:\n";
01253 printcpu();
01254 #endif
01255
01256
01257 for (i=0; i<=sizeL-winL; i++)
01258 {
01259 for (j=0; j<=sizeP-winP; j++)
01260 {
01261 real4 p = NORMS2(i,j).real() * NORMS2(i,j).imag();
01262 Result(i2+i,j2+j) = (p>0.0) ? CINT2(i,j)/sqrt(p) : 0.0;
01263 }
01264 }
01265
01266 #ifdef __DEBUGMAT2
01267 cout << "time for coherence.:\n";
01268 printcpu();
01269 #endif
01270
01271
01272 j2 += (SIZEBLOCKP - winP + 1);
01273 win.pixlo = win.pixhi + 1 - leadingzeros - trailingzeros;
01274 win.pixhi = win.pixlo + sizeP - 1;
01275 if (win.pixhi > CINT.pixels() - 1)
01276 {
01277 cout << "BERT: win.pixhi > CINT.pixels() - 1\n";
01278 if (win.pixlo+leadingzeros > Result.pixels()-trailingzeros-1)
01279 newline = true;
01280 if (newline)
01281 {
01282 cout << "BERT: newline==true\n";
01283 if (nonewline)
01284 break;
01285 newline = false;
01286 sizeP = SIZEBLOCKP;
01287 win.linelo = win.linehi + 1 - (winL-1)/2;
01288 win.linehi = win.linelo + sizeL - 1;
01289 win.pixlo = 0;
01290 win.pixhi = win.pixlo + sizeP - 1;
01291 i2 += (SIZEBLOCKL - winL + 1);
01292 j2 = leadingzeros;
01293 if (win.linehi > CINT.lines() - 1)
01294 {
01295 nonewline = true;
01296 win.linehi = CINT.lines() - 1;
01297 sizeL = win.linehi - win.linelo + 1;
01298 }
01299 }
01300
01301 else
01302 {
01303 cout << "BERT: newline==false\n";
01304 newline = true;
01305 win.pixhi = CINT.pixels() - 1;
01306 sizeP = win.pixhi - win.pixlo + 1;
01307 }
01308
01309 win2.linehi = sizeL - 1;
01310 win2.pixhi = sizeP - 1;
01311 }
01312 }
01313 #undef FFTS
01314 #endif
01315
01316
01317
01318
01319 return Result;
01320 }
01321
01322
01323
01324
01325
01326
01327
01328
01329
01330
01331 matrix<real4> coherence2(
01332 const matrix<complr4> &CINT,
01333 const matrix<complr4> &NORMS,
01334 uint winL,
01335 uint winP)
01336 {
01337 #ifdef __DEBUGMAT2
01338 matDEBUG.print("coherence2");
01339 if (!winL>=winP)
01340 matDEBUG.print("coherence: estimator window size L<P not very efficiently programmed.");
01341 #endif
01342 #ifdef __DEBUGMAT1
01343 if (CINT.lines() != NORMS.lines() ||
01344 CINT.pixels() != NORMS.pixels() )
01345 matERROR.print("coherence2::not same dimensions.");
01346 #endif
01347
01348 matrix<real4> Result(CINT.lines()-winL+1,CINT.pixels());
01349 int32 leadingzeros = (winP-1)/2;
01350 int32 trailingzeros = (winP)/2;
01351
01352 register int32 i,j,k,l;
01353 register complr4 sum;
01354 register complr4 power;
01355 for (j=leadingzeros; j<Result.pixels()-trailingzeros; j++)
01356 {
01357 sum = complr4(0.);
01358 power = complr4(0.);
01359
01360
01361 for (k=0; k<winL; k++)
01362 {
01363 for (l=j-leadingzeros; l<j-leadingzeros+winP; l++)
01364 {
01365 sum += CINT(k,l);
01366 power += NORMS(k,l);
01367 }
01368 }
01369
01370 real4 p = power.real() * power.imag();
01371 Result(0,j) = (p>0.0) ? sqrt(norm(sum)/p) : 0.0;
01372
01373
01374 for (i=0; i<Result.lines()-1; i++)
01375 {
01376 for (l=j-leadingzeros; l<j-leadingzeros+winP; l++)
01377 {
01378 sum += (CINT(i+winL,l) - CINT(i,l));
01379 power += (NORMS(i+winL,l) - NORMS(i,l));
01380 }
01381 real4 p = power.real() * power.imag();
01382 Result(i+1,j) = (p>0.0) ? sqrt(norm(sum)/p) : 0.0;
01383 }
01384 }
01385 return Result;
01386 }
01387
01388
01389
01390
01391
01392
01393
01394
01395
01396
01397
01398
01399 int32 mycomp2 (const void *x1, const void *x2)
01400 {
01401
01402
01403
01404 int ret = 0;
01405 if (*(float*)x1 < *(float*)x2 ) ret = -1;
01406 else if (*(float*)x1 > *(float*)x2 ) ret = 1;
01407 else if (*((float*)x1+1) < *((float*)x2+1) ) ret = -1;
01408 else if (*((float*)x1+1) > *((float*)x2+1) ) ret = 1;
01409
01410 return ret;
01411 }
01412
01413 void mysort2 (matrix<real4> &A)
01414 {
01415 #ifdef __DEBUGMAT2
01416 matDEBUG.print("mysort2 (real4)");
01417 #endif
01418 #ifdef __DEBUGMAT1
01419 if (A.pixels() < 2)
01420 matERROR.print("mysort sorts only min. 2 cols.");
01421 #endif
01422 qsort(&A[0][0],A.lines(),A.pixels()*sizeof(real4),mycomp2);
01423 }
01424
01425 void mysort2 (matrix<int32> &A)
01426 {
01427 #ifdef __DEBUGMAT2
01428 matDEBUG.print("mysort2 (int)");
01429 #endif
01430 #ifdef __DEBUGMAT1
01431 if (A.pixels() < 2)
01432 matERROR.print("mysort sorts only min. 2 cols.");
01433 #endif
01434 qsort(&A[0][0],A.lines(),A.pixels()*sizeof(int32),mycomp2);
01435 }
01436
01437
01438
01439
01440
01441
01442
01443
01444
01445
01446
01447
01448
01449
01450
01451
01452
01453
01454
01455
01456
01457
01458
01459
01460
01461
01462
01463
01464
01465
01466
01467
01468
01469
01470
01471
01472
01473
01474
01475
01476
01477
01478
01479
01480
01481
01482
01483
01484
01485
01486
01487
01488
01489
01490
01491
01492
01493
01494
01495
01496
01497
01498
01499
01500
01501
01502
01503
01504 #ifndef __USE_VECLIB_LIBRARY__
01505 #ifndef __USE_FFTW_LIBRARY__
01506 #define SWAP(a,b) bert=(a);(a)=(b);(b)=bert
01507 #endif
01508 #endif
01509
01510 void four1(
01511 complr4 data[],
01512 int32 fftlength,
01513 int32 isign)
01514 {
01515
01516
01517
01518
01519 #ifdef __DEBUGMAT1
01520 if (!myispower2(fftlength))
01521 matERROR.print("four1: length fft must be power of 2");
01522 if (abs(isign) != 1)
01523 matERROR.print("four1: isign should be (-)1.");
01524 #endif
01525
01526
01527 #ifdef __USE_FFTW_LIBRARY__
01528
01529
01530
01531
01532
01533
01534
01535
01536 #ifdef __NO_GURU_YET__
01537 fftwf_plan fftwf_tmp_plan = (isign==1) ?
01538 fftwf_plan_dft_1d(fftlength,
01539 reinterpret_cast<fftwf_complex*>(data),
01540 reinterpret_cast<fftwf_complex*>(data),
01541 FFTW_FORWARD, FFTW_ESTIMATE | FFTW_DESTROY_INPUT) :
01542 fftwf_plan_dft_1d(fftlength, r
01543 reinterpret_cast<fftwf_complex*>(data),
01544 reinterpret_cast<fftwf_complex*>(data),
01545 FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_DESTROY_INPUT) :
01546 fftwf_execute(fftwf_tmp_plan);
01547 if (isign==-1)
01548 for (int32 i=0; i<fftlength; ++i)
01549 #ifdef __GNUC__
01550 data[i] /= fftlength;
01551 #else
01552 data[i] /= complr4(fftlength);
01553 #endif
01554 fftwf_destroy_plan(fftwf_tmp_plan);
01555
01556
01557
01558 #else
01559 if (isign == 1)
01560 {
01561 static int32 last_length_fwd = 0;
01562 static fftwf_plan fftwf_fwd_plan;
01563 if (fftlength != last_length_fwd)
01564 {
01565 #ifdef __DEBUGMAT2
01566 matDEBUG.print("BK: FFTW: creating guru plan...");
01567 matDEBUG << "BK: FFTW: fftlength: " << fftlength
01568 << "; last_length_fwd: " << last_length_fwd;
01569 matDEBUG.print();
01570 #endif
01571 last_length_fwd = fftlength;
01572 fftwf_fwd_plan = fftwf_plan_dft_1d(fftlength,
01573 reinterpret_cast<fftwf_complex*>(data),
01574 reinterpret_cast<fftwf_complex*>(data),
01575 FFTW_FORWARD,
01576 FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
01577 #ifdef __DEBUGMAT2
01578 fftwf_print_plan(fftwf_fwd_plan);
01579 matDEBUG.print(" ");
01580 #endif
01581 }
01582
01583 #ifdef __DEBUGMAT2
01584 matDEBUG << "BK: FFTW: executing plan..." << "last_length_fwd: " << last_length_fwd;
01585 matDEBUG.print();
01586 #endif
01587 fftwf_execute_dft(fftwf_fwd_plan,
01588 reinterpret_cast<fftwf_complex*>(data),
01589 reinterpret_cast<fftwf_complex*>(data));
01590 }
01591
01592
01593 else
01594 {
01595 static int32 last_length_inv = 0;
01596 static fftwf_plan fftwf_inv_plan;
01597 if (fftlength != last_length_inv)
01598 {
01599 #ifdef __DEBUGMAT2
01600 matDEBUG.print("BK: FFTW: creating guru plan...");
01601 matDEBUG << "BK: FFTW: fftlength: " << fftlength
01602 << "; last_length_inv: " << last_length_inv;
01603 matDEBUG.print();
01604 #endif
01605 last_length_inv = fftlength;
01606 fftwf_inv_plan = fftwf_plan_dft_1d(fftlength,
01607 reinterpret_cast<fftwf_complex*>(data),
01608 reinterpret_cast<fftwf_complex*>(data),
01609 FFTW_BACKWARD,
01610 FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
01611 #ifdef __DEBUGMAT2
01612 fftwf_print_plan(fftwf_inv_plan);
01613 matDEBUG.print(" ");
01614 #endif
01615 }
01616
01617 #ifdef __DEBUGMAT2
01618 matDEBUG << "BK: FFTW: executing plan..." << "last_length_inv: " << last_length_inv;
01619 matDEBUG.print();
01620 #endif
01621 fftwf_execute_dft(fftwf_inv_plan,
01622 reinterpret_cast<fftwf_complex*>(data),
01623 reinterpret_cast<fftwf_complex*>(data));
01624
01625 for (int32 i=0; i<fftlength; ++i)
01626 #ifdef __GNUC__
01627 data[i] /= fftlength;
01628 #else
01629 data[i] /= complr4(fftlength);
01630 #endif
01631 }
01632 }
01633 #endif // guru interface for an unknown number of transforms
01634 #endif //fftw
01635
01636
01637
01638
01639
01640
01641
01642
01643
01644
01645 #ifndef __USE_FFTW_LIBRARY__
01646 #ifndef __USE_VECLIB_LIBRARY__
01647 register int32 i,j,n,mmax,m,istep;
01648 register complr4 bert;
01649 register double theta;
01650 register complr8 w,wp;
01651
01652 #define SAMEASVECLIB
01653 #ifdef SAMEASVECLIB
01654
01655
01656
01657
01658
01659 if (isign == -1)
01660 {for (i=1; i<fftlength/2; ++i)
01661 {SWAP(data[i],data[fftlength-i]);}}
01662 #endif
01663
01664 n = fftlength << 1;
01665 j = 1;
01666
01667
01668
01669
01670
01671 for (i=1; i<n; i+=2)
01672 {
01673 if (j > i)
01674 {
01675 SWAP(data[j/2],data[i/2]);
01676 }
01677 m = n >> 1;
01678 while (m>=2 && j>m)
01679 {
01680 j -= m;
01681 m >>= 1;
01682 }
01683 j += m;
01684 }
01685 mmax=2;
01686 while (n > mmax)
01687 {
01688 istep = mmax << 1;
01689 theta = isign*(6.28318530717959/mmax);
01690 wp = complr8(-2.*sqr(sin(.5*theta)),sin(theta));
01691 w = complr8(1.,0.);
01692 for (m=1; m<mmax; m+=2)
01693 {
01694 for (i=m; i<=n; i+=istep)
01695 {
01696 j = i+mmax;
01697 bert = data[j/2];
01698
01699 #ifdef __GNUC__
01700 bert *= complr4(real4(w.real()),real4(w.imag()));
01701 #else // assume aCC hp or so...
01702 bert *= w;
01703 #endif
01704 data[j/2] = data[i/2] - bert;
01705 data[i/2] += bert;
01706 }
01707 w += w * wp;
01708 }
01709 mmax = istep;
01710 }
01711
01712 #ifdef SAMEASVECLIB
01713
01714
01715
01716
01717 if (isign == 1)
01718 {for (i=1; i<fftlength/2; ++i)
01719 {SWAP(data[i],data[fftlength-i]);}}
01720 #endif
01721
01722
01723
01724 if (isign == -1)
01725 for (i=0; i<fftlength; ++i)
01726 #ifdef __GNUC__
01727 data[i] /= fftlength;
01728 #else
01729 data[i] /= complr4(fftlength);
01730 #endif
01731 #undef SWAP
01732 }
01733 #endif //__USE_VECLIB_LIBRARY__
01734 #endif //__USE_FFTW_LIBRARY__
01735
01736
01737
01738
01739 #ifndef __USE_FFTW_LIBRARY__
01740 #ifdef __USE_VECLIB_LIBRARY__
01741
01742
01743
01744
01745
01746
01747
01748
01749
01750
01751
01752 int32 ierr = 0;
01753 static matrix<real4> work;
01754 if (int32(2.5*fftlength) != work.size())
01755 {
01756 work.resize(1,int32(2.5*fftlength));
01757
01758 int32 iopt = -3;
01759 c1dfft(data,&fftlength,work[0],&iopt,&ierr);
01760 #ifdef __DEBUGMAT1
01761 switch (ierr)
01762 {
01763 case 0: matDEBUG.print("c1dfft ok."); break;
01764 case -1: matERROR.print("length < 0"); break;
01765 case -2: matERROR.print("l not of required form"); break;
01766 case -3: matERROR.print("invalid value of iopt"); break;
01767 case -4: matERROR.print("insufficient dynamical memory available in work"); break;
01768 default: matERROR.print("four1: unrecognized ierr.");
01769 }
01770 #endif
01771 }
01772
01773 c1dfft(data,&fftlength,work[0],&isign,&ierr);
01774 if (ierr != 0) {cerr << "veclib: c1dfft: ierr = " << ierr << endl; exit(-1);}
01775 }
01776 #endif // either internal or veclib complex 1d fft or fftw
01777 #endif // either internal or veclib complex 1d fft or fftw
01778
01779
01780
01781
01782
01783
01784
01785
01786
01787
01788
01789
01790
01791
01792
01793
01794
01795
01796
01797
01798
01799
01800
01801
01802
01803
01804
01805
01806 void fft(matrix<complr4> &A, int32 dimension)
01807 {
01808 register int32 i;
01809 const int32 iopt = 1;
01810 switch (dimension)
01811 {
01812 case 1:
01813 {
01814 #ifdef __DEBUGMAT2
01815 matDEBUG.print("1d fft over columns");
01816 #endif
01817 const int32 fftlength = A.lines();
01818 for (i=0; i<A.pixels(); ++i)
01819 {
01820
01821 matrix<complr4> VECTOR = A.getcolumn(i);
01822 four1(&VECTOR[0][0],fftlength,iopt);
01823 A.setcolumn(i,VECTOR);
01824 }
01825 break;
01826 }
01827 case 2:
01828 #ifdef __DEBUGMAT2
01829 matDEBUG.print("1d fft over rows");
01830 #endif
01831 {
01832 const int32 fftlength = A.pixels();
01833 for (i=0; i<A.lines(); ++i)
01834 four1(&A[i][0],fftlength,iopt);
01835 break;
01836 }
01837 default:
01838 matERROR.print("fft: dimension != {1,2}");
01839 }
01840 }
01841
01842
01843
01844
01845
01846
01847
01848
01849
01850
01851
01852
01853
01854 void ifft(matrix<complr4> &A, int32 dimension)
01855 {
01856 register int32 i;
01857 const int32 iopt = -1;
01858 switch (dimension)
01859 {
01860 case 1:
01861 {
01862 #ifdef __DEBUGMAT2
01863 matDEBUG.print("1d ifft over columns");
01864 #endif
01865 const int32 fftlength = A.lines();
01866 for (i=0; i<A.pixels(); ++i)
01867 {
01868 matrix<complr4> VECTOR = A.getcolumn(i);
01869 four1(&VECTOR[0][0],fftlength,iopt);
01870 A.setcolumn(i,VECTOR);
01871 }
01872 break;
01873 }
01874 case 2:
01875 #ifdef __DEBUGMAT2
01876 matDEBUG.print("1d ifft over rows");
01877 #endif
01878 {
01879 const int32 fftlength = A.pixels();
01880 for (i=0; i<A.lines(); ++i)
01881 four1(&A[i][0],fftlength,iopt);
01882 break;
01883 }
01884 default:
01885 matERROR.print("fft: dimension != {1,2}");
01886 }
01887 }
01888
01889
01890
01891
01892
01893
01894
01895
01896 void fft2d(
01897 matrix<complr4> &A)
01898 {
01899 #ifndef __USE_FFTW_LIBRARY__ // prefer FFTW
01900 #ifdef __USE_VECLIB_LIBRARY__ // use VECLIB
01901 #ifdef __DEBUGMAT2
01902 matDEBUG.print("fft2d: veclib");
01903 #endif
01904 int32 l2 = A.lines();
01905 int32 l1 = A.pixels();
01906 int32 ldz= l1;
01907 int32 iopt = 1;
01908 int32 ierr;
01909 c2dfft(A[0],&l1,&l2,&ldz,&iopt,&ierr);
01910 #ifdef __DEBUGMAT1
01911 switch (ierr)
01912 {
01913 case 0: matDEBUG.print("fft2d: ok"); break;
01914 case -1: matERROR.print("fft2d: l1 not of the required form"); break;
01915 case -2: matERROR.print("fft2d: l2 not of the required form"); break;
01916 case -3: matERROR.print("fft2d: ldz not of the required form"); break;
01917 case -4: matERROR.print("fft2d: probable error in ldz or dimension of z"); break;
01918 default: matERROR.print("fft2d: unrecognized error.");
01919 }
01920 #endif
01921 }
01922 #endif
01923 #endif
01924
01925
01926
01927 #ifdef __USE_FFTW_LIBRARY__ // prefer FFTW
01928
01929
01930
01931
01932
01933
01934
01935
01936 #ifdef __DEBUGMAT2
01937 matDEBUG.print("fft2d: fftw");
01938 #endif
01939 int32 nx = A.pixels();
01940 int32 ny = A.lines();
01941 fftwf_plan fftwf_tmp_plan;
01942 #ifdef __DEBUGMAT2
01943 matDEBUG.print("BK: FFTW: creating tmp fwd 2d plan.");
01944 #endif
01945 fftwf_tmp_plan = fftwf_plan_dft_2d(
01946 nx, ny,
01947 reinterpret_cast<fftwf_complex*>(A[0]),
01948 reinterpret_cast<fftwf_complex*>(A[0]),
01949 FFTW_FORWARD,
01950 FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
01951 #ifdef __DEBUGMAT2
01952 fftwf_print_plan(fftwf_tmp_plan);
01953 matDEBUG.print(" ");
01954 #endif
01955
01956 #ifdef __DEBUGMAT2
01957 matDEBUG.print("BK: FFTW: executing tmp fwd 2d plan.");
01958 #endif
01959 fftwf_execute(fftwf_tmp_plan);
01960 #ifdef __DEBUGMAT2
01961 matDEBUG.print("BK: FFTW: destroying tmp fwd 2d plan.");
01962 #endif
01963 fftwf_destroy_plan(fftwf_tmp_plan);
01964 }
01965 #endif // select FFTW
01966
01967
01968
01969 #ifndef __USE_FFTW_LIBRARY__
01970 #ifndef __USE_VECLIB_LIBRARY__
01971 #ifdef __DEBUGMAT2
01972 matDEBUG.print("fft2d: internal (using multiple 1d ffts)");
01973 #endif
01974 fft(A,2);
01975 fft(A,1);
01976 }
01977 #endif // select internal
01978 #endif // select internal
01979
01980
01981
01982
01983
01984
01985
01986
01987
01988 void ifft2d(
01989 matrix<complr4> &A)
01990 {
01991 #ifndef __USE_FFTW_LIBRARY__ // prefer fftw
01992 #ifdef __USE_VECLIB_LIBRARY__ // but use veclib before internal implementation.
01993 #ifdef __DEBUGMAT2
01994 matDEBUG.print("ifft2d: veclib");
01995 #endif
01996 int32 l2 = A.lines();
01997 int32 l1 = A.pixels();
01998 int32 ldz= l1;
01999 int32 iopt = -1;
02000 int32 ierr;
02001 c2dfft(A[0],&l1,&l2,&ldz,&iopt,&ierr);
02002 #ifdef __DEBUGMAT1
02003 switch (ierr)
02004 {
02005 case 0: matDEBUG.print("ifft2d: ok"); break;
02006 case -1: matERROR.print("ifft2d: l1 not of the required form"); break;
02007 case -2: matERROR.print("ifft2d: l2 not of the required form"); break;
02008 case -3: matERROR.print("ifft2d: ldz not of the required form"); break;
02009 case -4: matERROR.print("ifft2d: probable error in ldz or dimension of z"); break;
02010 default: matERROR.print("ifft2d: unrecognized error.");
02011 }
02012 #endif
02013 }
02014 #endif
02015 #endif
02016
02017
02018
02019 #ifdef __USE_FFTW_LIBRARY__ // prefer fftw
02020
02021
02022
02023
02024
02025
02026
02027 #ifdef __DEBUGMAT2
02028 matDEBUG.print("ifft2d: fftw");
02029 #endif
02030 int32 nx = A.pixels();
02031 int32 ny = A.lines();
02032 fftwf_plan fftwf_tmp_plan;
02033 #ifdef __DEBUGMAT2
02034 matDEBUG.print("BK: FFTW: creating tmp inv 2d plan.");
02035 #endif
02036 fftwf_tmp_plan = fftwf_plan_dft_2d(
02037 nx, ny,
02038 reinterpret_cast<fftwf_complex*>(A[0]),
02039 reinterpret_cast<fftwf_complex*>(A[0]),
02040 FFTW_BACKWARD,
02041 FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
02042 #ifdef __DEBUGMAT2
02043 fftwf_print_plan(fftwf_tmp_plan);
02044 matDEBUG.print(" ");
02045 #endif
02046
02047 #ifdef __DEBUGMAT2
02048 matDEBUG.print("BK: FFTW: executing tmp inv 2d plan.");
02049 #endif
02050 fftwf_execute(fftwf_tmp_plan);
02051 #ifdef __DEBUGMAT2
02052 matDEBUG.print("BK: FFTW: destroying tmp inv 2d plan.");
02053 #endif
02054 fftwf_destroy_plan(fftwf_tmp_plan);
02055
02056 #ifdef __DEBUGMAT2
02057 matDEBUG.print("BK: FFTW: normalizing inverse transform.");
02058 #endif
02059 A /= complr4(nx*ny);
02060
02061
02062
02063
02064
02065
02066
02067 }
02068 #endif
02069
02070
02071
02072
02073 #ifndef __USE_FFTW_LIBRARY__
02074 #ifndef __USE_VECLIB_LIBRARY__
02075 #ifdef __DEBUGMAT2
02076 matDEBUG.print("ifft2d: internal");
02077 #endif
02078 ifft(A,2);
02079 ifft(A,1);
02080 }
02081 #endif
02082 #endif
02083
02084
02085
02086
02087
02088
02089
02090
02091
02092
02093
02094
02095 matrix<complr4> oversample(
02096 matrix<complr4> A,
02097 uint factorrow,
02098 uint factorcol)
02099 {
02100 #ifdef __DEBUGMAT2
02101 matDEBUG.print("oversample");
02102 #endif
02103 const uint l = A.lines();
02104 const uint p = A.pixels();
02105 const uint halfl = l/2;
02106 const uint halfp = p/2;
02107 const uint L2 = factorrow*l;
02108 const uint P2 = factorcol*p;
02109
02110 #ifdef __DEBUGMAT1
02111 if (A.isvector())
02112 matERROR.print("OVERSAMPLE: only 2d matrices.");
02113 if (!myispower2(l) && factorrow != 1)
02114 matERROR.print("OVERSAMPLE: numlines != 2^n.");
02115 if (!myispower2(p) && factorcol != 1)
02116 matERROR.print("OVERSAMPLE: numcols != 2^n.");
02117 #endif
02118
02119 register int32 i,j;
02120 const real4 half = 0.5;
02121 matrix<complr4> Res(L2,P2);
02122 if (factorrow==1)
02123 {
02124 fft(A,2);
02125 for (i=0; i<l; ++i)
02126 #ifdef __GNUC__
02127 A(i,halfp) *= half;
02128 #else
02129 A(i,halfp) *= complr4(half);
02130 #endif
02131 const window winA1(0, l-1, 0, halfp);
02132 const window winA2(0, l-1, halfp, p-1);
02133 const window winR2(0, l-1, P2-halfp, P2-1);
02134 Res.setdata(winA1,A,winA1);
02135 Res.setdata(winR2,A,winA2);
02136 ifft(Res,2);
02137 }
02138 else if (factorcol==1)
02139 {
02140 fft(A,1);
02141 for (i=0; i<p; ++i)
02142 #ifdef __GNUC__
02143 A(halfl,i) *= half;
02144 #else
02145 A(halfl,i) *= complr4(half);
02146 #endif
02147 const window winA1(0, halfl, 0, p-1);
02148 const window winA2(halfl, l-1, 0, p-1);
02149 const window winR2(L2-halfl, L2-1, 0, p-1);
02150 Res.setdata(winA1,A,winA1);
02151 Res.setdata(winR2,A,winA2);
02152 ifft(Res,1);
02153 }
02154 else
02155 {
02156 fft2d(A);
02157 for (i=0; i<l; ++i)
02158 #ifdef __GNUC__
02159 A(i,halfp) *= half;
02160 #else
02161 A(i,halfp) *= complr4(half);
02162 #endif
02163 for (i=0; i<p; ++i)
02164 #ifdef __GNUC__
02165 A(halfl,i) *= half;
02166 #else
02167 A(halfl,i) *= complr4(half);
02168 #endif
02169 const window winA1(0, halfl, 0, halfp);
02170 const window winA2(0, halfl, halfp, p-1);
02171 const window winA3(halfl, l-1, 0, halfp);
02172 const window winA4(halfl, l-1, halfp, p-1);
02173 const window winR2(0, halfl, P2-halfp, P2-1);
02174 const window winR3(L2-halfl, L2-1, 0, halfp);
02175 const window winR4(L2-halfl, L2-1, P2-halfp, P2-1);
02176 Res.setdata(winA1,A,winA1);
02177 Res.setdata(winR2,A,winA2);
02178 Res.setdata(winR3,A,winA3);
02179 Res.setdata(winR4,A,winA4);
02180 ifft2d(Res);
02181 }
02182
02183 Res *= real4(factorrow*factorcol);
02184 return Res;
02185
02186 }
02187
02188
02189
02190
02191
02192
02193
02194 matrix<real4> oversample(
02195 const matrix<real4> &A,
02196 uint factorrow,
02197 uint factorcol)
02198 {
02199
02200 matrix<complr4> TMP = oversample(mat2cr4(A),factorrow,factorcol);
02201 return real(TMP);
02202 }
02203
02204
02205
02206
02207
02208
02209
02210
02211 void dotmult(
02212 complr4 *pntR,
02213 const matrix<real4> &B,
02214 int32 stride)
02215 {
02216 #ifdef __DEBUGMAT2
02217 matDEBUG.print("dotmult, no range checking");
02218 #endif
02219
02220
02221
02222
02223 real4 *pntB = B[0];
02224 for (register int32 i=0; i<B.size(); i++)
02225 {
02226 #ifdef __GNUC__
02227 (*pntR) *= (*pntB++);
02228 #else
02229 (*pntR) *= complr4(*pntB++);
02230 #endif
02231 pntR += stride;
02232 }
02233 }
02234
02235
02236
02237
02238
02239 #ifdef __USE_VECLIB_LIBRARY__
02240
02241
02242
02243
02244 matrix<real4> operator * (const matrix<real4>& A, const matrix<real4>& B)
02245 {
02246 #ifdef __DEBUGMAT2
02247 matDEBUG.print("matrices * (veclib)");
02248 #endif
02249 #ifdef __DEBUGMAT1
02250 if (A.pixels() != B.lines())
02251 matERROR.print("matrix::operator *: multiplication not possible");
02252 if (!A.size())
02253 matERROR.print("matrixbaseclass: operator * with empty matrices.");
02254 #endif
02255 matrix<real4> Result(A.lines(),B.pixels());
02256 int32 n = B.pixels();
02257 int32 m = A.lines();
02258 int32 k = A.pixels();
02259 real4 r4alpha = 1.;
02260 real4 r4beta = 0.;
02261 sgemm("N","N",&n,&m,&k,&r4alpha,B[0],&n,A[0],&k,&r4beta,Result[0],&n,1,1);
02262 return Result;
02263 }
02264 #endif
02265
02266
02267 #ifdef __USE_VECLIB_LIBRARY__
02268
02269
02270
02271
02272 matrix<real8> operator * (const matrix<real8>& A, const matrix<real8>& B)
02273 {
02274 #ifdef __DEBUGMAT2
02275 matDEBUG.print("matrices * (veclib)");
02276 #endif
02277 #ifdef __DEBUGMAT1
02278 if (A.pixels() != B.lines())
02279 matERROR.print("matrix::operator *: multiplication not possible");
02280 if (!A.size())
02281 matERROR.print("matrixbaseclass: operator * with empty matrices.");
02282 #endif
02283 matrix<real8> Result(A.lines(),B.pixels());
02284 int32 n = B.pixels();
02285 int32 m = A.lines();
02286 int32 k = A.pixels();
02287 real8 r8alpha =1.;
02288 real8 r8beta = 0.;
02289 dgemm("N","N",&n,&m,&k,&r8alpha,B[0],&n,A[0],&k,&r8beta,Result[0],&n,1,1);
02290 return Result;
02291 }
02292 #endif
02293
02294
02295 #ifdef __USE_VECLIB_LIBRARY__
02296
02297
02298
02299
02300 matrix<complr4> operator * (const matrix<complr4>& A, const matrix<complr4>& B)
02301 {
02302 #ifdef __DEBUGMAT2
02303 matDEBUG.print("matrices * (veclib)");
02304 #endif
02305 #ifdef __DEBUGMAT1
02306 if (A.pixels() != B.lines())
02307 matERROR.print("matrix::operator *: multiplication not possible");
02308 if (!A.size())
02309 matERROR.print("matrixbaseclass: operator * with empty matrices.");
02310 #endif
02311 matrix<complr4> Result(A.lines(),B.pixels());
02312 int32 n = B.pixels();
02313 int32 m = A.lines();
02314 int32 k = A.pixels();
02315 complr4 c4alpha(1.,0.);
02316 complr4 c4beta(0.,0.);
02317 cgemm("N","N",&n,&m,&k,&c4alpha,B[0],&n,A[0],&k,&c4beta,Result[0],&n,1,1);
02318 return Result;
02319 }
02320 #endif
02321
02322
02323 #ifdef __USE_VECLIB_LIBRARY__
02324
02325
02326
02327
02328 matrix<complr8> operator * (const matrix<complr8>& A, const matrix<complr8>& B)
02329 {
02330 #ifdef __DEBUGMAT2
02331 matDEBUG.print("matrices * (veclib)");
02332 #endif
02333 #ifdef __DEBUGMAT1
02334 if (A.pixels() != B.lines())
02335 matERROR.print("matrix::operator *: multiplication not possible");
02336 if (!A.size())
02337 matERROR.print("matrixbaseclass: operator * with empty matrices.");
02338 #endif
02339 matrix<complr8> Result(A.lines(),B.pixels());
02340 int32 n = B.pixels();
02341 int32 m = A.lines();
02342 int32 k = A.pixels();
02343 complr8 c8alpha = 1.;
02344 complr8 c8beta = 0.;
02345 zgemm("N","N",&n,&m,&k,&c8alpha,B[0],&n,A[0],&k,&c8beta,Result[0],&n,1,1);
02346 return Result;
02347 }
02348 #endif
02349
02350
02351
02352 #ifdef __USE_VECLIB_LIBRARY__
02353
02354
02355
02356
02357 matrix<real4> matTxmat(const matrix<real4> &A, const matrix<real4> &B)
02358 {
02359 #ifdef __DEBUGMAT2
02360 matDEBUG.print("matTxmat, veclib::sgemm");
02361 #endif
02362 #ifdef __DEBUGMAT1
02363 if (A.lines() != B.lines())
02364 matERROR.print("matTxmat: size A,B: input is A,B; computed is trans(A)*B.");
02365 #endif
02366 matrix<real4> Result(A.pixels(),B.pixels());
02367 int32 m = A.pixels();
02368 int32 n = B.pixels();
02369 int32 k = B.lines();
02370 real4 r4alpha =1.;
02371 real4 r4beta = 0.;
02372 sgemm("N","T",&n,&m,&k,&r4alpha,B[0],&n,A[0],&m,&r4beta,Result[0],&n,1,1);
02373 return Result;
02374 }
02375 #endif
02376
02377
02378
02379 #ifdef __USE_VECLIB_LIBRARY__
02380
02381
02382
02383
02384 matrix<real8> matTxmat(const matrix<real8> &A, const matrix<real8> &B)
02385 {
02386 #ifdef __DEBUGMAT2
02387 matDEBUG.print("matTxmat, veclib::dgemm");
02388 #endif
02389 #ifdef __DEBUGMAT1
02390 if (A.lines() != B.lines())
02391 matERROR.print("matTxmat: size A,B: input is A,B; computed is trans(A)*B.");
02392 #endif
02393 matrix<real8> Result(A.pixels(),B.pixels());
02394 int32 m = A.pixels();
02395 int32 n = B.pixels();
02396 int32 k = B.lines();
02397 real8 r8alpha = 1.;
02398 real8 r8beta = 0.;
02399 dgemm("N","T",&n,&m,&k,&r8alpha,B[0],&n,A[0],&m,&r8beta,Result[0],&n,1,1);
02400 return Result;
02401 }
02402 #endif
02403
02404
02405
02406 #ifdef __USE_VECLIB_LIBRARY__
02407
02408
02409
02410
02411 matrix<complr4> matTxmat(const matrix<complr4> &A, const matrix<complr4> &B)
02412 {
02413 #ifdef __DEBUGMAT2
02414 matDEBUG.print("matTxmat, veclib::cgemm");
02415 #endif
02416 #ifdef __DEBUGMAT1
02417 if (A.lines() != B.lines())
02418 matERROR.print("matTxmat: size A,B: input is A,B; computed is trans(A)*B.");
02419 #endif
02420 matrix<complr4> Result(A.pixels(),B.pixels());
02421 int32 m = A.pixels();
02422 int32 n = B.pixels();
02423 int32 k = B.lines();
02424 complr4 c4alpha(1.,0.);
02425 complr4 c4beta(0.,0.);
02426 cgemm("N","T",&n,&m,&k,&c4alpha,B[0],&n,A[0],&m,&c4beta,Result[0],&n,1,1);
02427 return Result;
02428 }
02429 #endif
02430
02431
02432
02433 #ifdef __USE_VECLIB_LIBRARY__
02434
02435
02436
02437
02438 matrix<complr8> matTxmat(const matrix<complr8> &A, const matrix<complr8> &B)
02439 {
02440 #ifdef __DEBUGMAT2
02441 matDEBUG.print("matTxmat, veclib::zgemm");
02442 #endif
02443 #ifdef __DEBUGMAT1
02444 if (A.lines() != B.lines())
02445 matERROR.print("matTxmat: size A,B: input is A,B; computed is trans(A)*B.");
02446 #endif
02447 matrix<complr8> Result(A.pixels(),B.pixels());
02448 int32 m = A.pixels();
02449 int32 n = B.pixels();
02450 int32 k = B.lines();
02451 complr8 c8alpha = 1.;
02452 complr8 c8beta = 0.;
02453 zgemm("N","T",&n,&m,&k,&c8alpha,B[0],&n,A[0],&m,&c8beta,Result[0],&n,1,1);
02454 return Result;
02455 }
02456 #endif
02457
02458
02459
02460 #ifdef __USE_VECLIB_LIBRARY__
02461
02462
02463
02464
02465 matrix<real4> matxmatT(const matrix<real4> &A, const matrix<real4> &B)
02466 {
02467 #ifdef __DEBUGMAT2
02468 matDEBUG.print("matxmatT, veclib::sgemm");
02469 #endif
02470 #ifdef __DEBUGMAT1
02471 if (A.pixels() != B.pixels())
02472 matERROR.print("matxmatT: size A,B: input is A,B; computed is A*trans(B).");
02473 #endif
02474 matrix<real4> Result(A.lines(),B.lines());
02475 int32 m = A.lines();
02476 int32 n = B.lines();
02477 int32 k = B.pixels();
02478 real4 r4alpha =1.;
02479 real4 r4beta = 0.;
02480 sgemm("T","N",&n,&m,&k,&r4alpha,B[0],&k,A[0],&k,&r4beta,Result[0],&n,1,1);
02481 return Result;
02482 }
02483 #endif
02484
02485
02486
02487 #ifdef __USE_VECLIB_LIBRARY__
02488
02489
02490
02491
02492 matrix<real8> matxmatT(const matrix<real8> &A, const matrix<real8> &B)
02493 {
02494 #ifdef __DEBUGMAT2
02495 matDEBUG.print("matxmatT, veclib::dgemm");
02496 #endif
02497 #ifdef __DEBUGMAT1
02498 if (A.pixels() != B.pixels())
02499 matERROR.print("matxmatT: size A,B: input is A,B; computed is A*trans(B).");
02500 #endif
02501 matrix<real8> Result(A.lines(),B.lines());
02502 int32 m = A.lines();
02503 int32 n = B.lines();
02504 int32 k = B.pixels();
02505 real8 r8alpha =1.;
02506 real8 r8beta = 0.;
02507 dgemm("T","N",&n,&m,&k,&r8alpha,B[0],&k,A[0],&k,&r8beta,Result[0],&n,1,1);
02508 return Result;
02509 }
02510 #endif
02511
02512
02513
02514 #ifdef __USE_VECLIB_LIBRARY__
02515
02516
02517
02518
02519 matrix<complr4> matxmatT(const matrix<complr4> &A, const matrix<complr4> &B)
02520 {
02521 #ifdef __DEBUGMAT2
02522 matDEBUG.print("matxmatT, veclib::cgemm");
02523 #endif
02524 #ifdef __DEBUGMAT1
02525 if (A.pixels() != B.pixels())
02526 matERROR.print("matxmatT: size A,B: input is A,B; computed is A*trans(B).");
02527 #endif
02528 matrix<complr4> Result(A.lines(),B.lines());
02529 int32 m = A.lines();
02530 int32 n = B.lines();
02531 int32 k = B.pixels();
02532 complr4 c4alpha(1.,0.);
02533 complr4 c4beta(0.,0.);
02534 cgemm("T","N",&n,&m,&k,&c4alpha,B[0],&k,A[0],&k,&c4beta,Result[0],&n,1,1);
02535 return Result;
02536 }
02537 #endif
02538
02539
02540
02541 #ifdef __USE_VECLIB_LIBRARY__
02542
02543
02544
02545
02546 matrix<complr8> matxmatT(const matrix<complr8> &A, const matrix<complr8> &B)
02547 {
02548 #ifdef __DEBUGMAT2
02549 matDEBUG.print("matxmatT, veclib::zgemm");
02550 #endif
02551 #ifdef __DEBUGMAT1
02552 if (A.pixels() != B.pixels())
02553 matERROR.print("matxmatT: size A,B: input is A,B; computed is A*trans(B).");
02554 #endif
02555 matrix<complr8> Result(A.lines(),B.lines());
02556 int32 m = A.lines();
02557 int32 n = B.lines();
02558 int32 k = B.pixels();
02559 complr8 c8alpha = 1.;
02560 complr8 c8beta = 0.;
02561 zgemm("T","N",&n,&m,&k,&c8alpha,B[0],&k,A[0],&k,&c8beta,Result[0],&n,1,1);
02562 return Result;
02563 }
02564 #endif
02565
02566
02567
02568
02569
02570
02571
02572 matrix<complr4>& matrix<complr4>::operator *= (const matrix<real4> &A)
02573 {
02574 #ifdef __DEBUGMAT2
02575 matDEBUG.print("*=");
02576 #endif
02577 #ifdef __DEBUGMAT1
02578 if (nrows != A.lines() || ncols != A.pixels())
02579 matERROR.print("matrixbaseclass: *= matrices must be same size.");
02580 #endif
02581 complr4 *pntmat = data[0];
02582 real4 *pntA = A[0];
02583 for (register int32 i=0; i<nsize; i++)
02584 (*pntmat++) *= (*pntA++);
02585 return *this;
02586 }
02587
02588
02589
02590
02591
02592
02593
02594 matrix<complr4>& matrix<complr4>::operator /= (const matrix<real4> &A)
02595 {
02596 #ifdef __DEBUGMAT2
02597 matDEBUG.print("/=");
02598 #endif
02599 #ifdef __DEBUGMAT1
02600 if (nrows != A.lines() || ncols != A.pixels())
02601 matERROR.print("matrixbaseclass: *= matrices must be same size.");
02602 #endif
02603 complr4 *pntmat = data[0];
02604 real4 *pntA = A[0];
02605 for (register int32 i=0; i<nsize; i++)
02606 (*pntmat++) /= (*pntA++);
02607 return *this;
02608 }
02609
02610
02611
02612
02613
02614
02615
02616 matrix<complr4>& matrix<complr4>::operator += (const matrix<real4> &A)
02617 {
02618 #ifdef __DEBUGMAT2
02619 matDEBUG.print("+=");
02620 #endif
02621 #ifdef __DEBUGMAT1
02622 if (nrows != A.lines() || ncols != A.pixels())
02623 matERROR.print("matrixbaseclass: += matrices must be same size.");
02624 #endif
02625 complr4 *pntmat = data[0];
02626 real4 *pntA = A[0];
02627 for (register int32 i=0; i<nsize; i++)
02628 (*pntmat++) += (*pntA++);
02629 return *this;
02630 }
02631
02632
02633
02634
02635
02636
02637
02638 matrix<complr4>& matrix<complr4>::operator -= (const matrix<real4> &A)
02639 {
02640 #ifdef __DEBUGMAT2
02641 matDEBUG.print("-=");
02642 #endif
02643 #ifdef __DEBUGMAT1
02644 if (nrows != A.lines() || ncols != A.pixels())
02645 matERROR.print("matrixbaseclass: -= matrices must be same size.");
02646 #endif
02647 complr4 *pntmat = data[0];
02648 real4 *pntA = A[0];
02649 for (register int32 i=0; i<nsize; i++)
02650 (*pntmat++) -= (*pntA++);
02651 return *this;
02652 }
02653
02654
02655
02656
02657
02658
02659
02660 matrix<complr4> diagxmat (const matrix<real4> &diag, const matrix<complr4> &B)
02661 {
02662 #ifdef __DEBUGMAT2
02663 matDEBUG.print("diagxmat: R4*CR4");
02664 #endif
02665 #ifdef __DEBUGMAT1
02666 if (min(diag.lines(),diag.pixels()) != 1)
02667 matERROR.print("diagxmat: sizes A,B: diag is vector.");
02668 if (diag.size() != B.lines())
02669 matERROR.print("diagxmat: sizes A,B: input is vector, matrix.");
02670 #endif
02671
02672 matrix<complr4> Result=B;
02673 if (diag.lines() != 1)
02674 {
02675 for (register int32 i=0; i<Result.lines(); i++)
02676 for (register int32 j=0; j<Result.pixels(); j++)
02677 Result(i,j) *= diag(i,0);
02678
02679 }
02680 else
02681 {
02682 for (register int32 i=0; i<Result.lines(); i++)
02683 for (register int32 j=0; j<Result.pixels(); j++)
02684 Result(i,j) *= diag(0,i);
02685
02686 }
02687 return Result;
02688 }
02689
02690
02691
02692
02693
02694
02695
02696 matrix<real4> cos(
02697 const matrix<real4> &A)
02698 {
02699 #ifdef __DEBUGMAT2
02700 matDEBUG.print("cos");
02701 #endif
02702 matrix<real4> Result(A.lines(),A.pixels());
02703 real4 *pntR = Result[0];
02704 real4 *pntA = A[0];
02705
02706 for (register int32 i=0; i<A.size(); i++)
02707 (*pntR++) = cos(*pntA++);
02708 return Result;
02709 }
02710
02711
02712
02713
02714
02715
02716
02717 matrix<real8> cos(
02718 const matrix<real8> &A)
02719 {
02720 #ifdef __DEBUGMAT2
02721 matDEBUG.print("cos");
02722 #endif
02723 matrix<real8> Result(A.lines(),A.pixels());
02724 real8 *pntR = Result[0];
02725 real8 *pntA = A[0];
02726
02727 for (register int32 i=0; i<A.size(); i++)
02728 (*pntR++) = cos(*pntA++);
02729 return Result;
02730 }
02731
02732
02733
02734
02735
02736
02737
02738 matrix<real4> sin(
02739 const matrix<real4> &A)
02740 {
02741 #ifdef __DEBUGMAT2
02742 matDEBUG.print("sin");
02743 #endif
02744 matrix<real4> Result(A.lines(),A.pixels());
02745 real4 *pntR = Result[0];
02746 real4 *pntA = A[0];
02747
02748 for (register int32 i=0; i<A.size(); i++)
02749 (*pntR++) = sin(*pntA++);
02750 return Result;
02751 }
02752
02753
02754
02755
02756
02757
02758
02759 matrix<real8> sin(
02760 const matrix<real8> &A)
02761 {
02762 #ifdef __DEBUGMAT2
02763 matDEBUG.print("sin");
02764 #endif
02765 matrix<real8> Result(A.lines(),A.pixels());
02766 real8 *pntR = Result[0];
02767 real8 *pntA = A[0];
02768
02769 for (register int32 i=0; i<A.size(); i++)
02770 (*pntR++) = sin(*pntA++);
02771 return Result;
02772 }
02773
02774
02775
02776
02777
02778
02779
02780 matrix<complr4> mat2cr4(
02781 const matrix<real4> &A)
02782 {
02783 #ifdef __DEBUGMAT2
02784 matDEBUG.print("mat2cr4. (real4)");
02785 #endif
02786 matrix<complr4> Result(A.lines(),A.pixels());
02787 real4 *pntA = A[0];
02788 complr4 *pntR = Result[0];
02789 for (register int32 i=0; i<Result.size(); i++)
02790 (*pntR++) = complr4(*pntA++);
02791 return Result;
02792 }
02793
02794
02795
02796
02797
02798
02799
02800 matrix<complr4> mat2cr4(
02801 const matrix<real4>& A,
02802 const matrix<real4>& B)
02803 {
02804 #ifdef __DEBUGMAT1
02805 if (A.lines()!=B.lines() || A.lines()!=B.lines())
02806 matERROR.print("operator complr4, input matrices not same size");
02807 #endif
02808 #ifdef __DEBUGMAT2
02809 matDEBUG.print("operator complr4(r4,r4)");
02810 #endif
02811 matrix<complr4> Result(A.lines(),A.pixels());
02812 real4 *pntA = A[0];
02813 real4 *pntB = B[0];
02814 complr4 *pntR = Result[0];
02815 for (register int32 i=0; i<Result.size(); i++)
02816 (*pntR++) = complr4(*pntA++,*pntB++);
02817 return Result;
02818 }
02819
02820
02821
02822
02823
02824
02825
02826 matrix<complr4> mat2cr4(
02827 const matrix<real8>& A,
02828 const matrix<real8>& B)
02829 {
02830 #ifdef __DEBUGMAT1
02831 if (A.lines()!=B.lines() || A.lines()!=B.lines())
02832 matERROR.print("operator complr4, input matrices not same size");
02833 #endif
02834 #ifdef __DEBUGMAT2
02835 matDEBUG.print("operator complr4(r8,r8)");
02836 #endif
02837 matrix<complr4> Result(A.lines(),A.pixels());
02838 real8 *pntA = A[0];
02839 real8 *pntB = B[0];
02840 complr4 *pntR = Result[0];
02841 for (register int32 i=0; i<Result.size(); i++)
02842 (*pntR++) = complr4(*pntA++,*pntB++);
02843 return Result;
02844 }
02845
02846
02847
02848
02849
02850
02851
02852
02853
02854
02855
02856
02857
02858
02859
02860
02861
02862
02863
02864
02865
02866
02867
02868
02869
02870
02871
02872
02873
02874
02875
02876
02877
02878
02879
02880
02881
02882
02883
02884
02885
02886
02887
02888