Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Class Members | File Members

matrixspecs.cc

Go to the documentation of this file.
00001 /*
00002  * Copyright (c) 1999-2003 Bert Kampes
00003  * Copyright (c) 1999-2003 Delft University of Technology, The Netherlands
00004  *
00005  * This file is part of Doris, the Delft o-o radar interferometric software.
00006  *
00007  * Doris program is free software; you can redistribute it and/or modify
00008  * it under the terms of the GNU General Public License as published by
00009  * the Free Software Foundation; either version 2 of the License, or
00010  * (at your option) any later version.
00011  *
00012  * Doris is distributed in the hope that it will be useful,
00013  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00014  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015  * GNU General Public License for more details.
00016  *
00017  * You should have received a copy of the GNU General Public License
00018  * along with this program; if not, write to the Free Software
00019  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00020  *
00021  * Publications that contain results produced by the Doris software should
00022  * contain an acknowledgment. (For example: The interferometric processing
00023  * was performed using the freely available Doris software package developed
00024  * by the Delft Institute for Earth-Oriented Space Research (DEOS), Delft
00025  * University of Technology, or include a reference to: Bert Kampes and
00026  * Stefania Usai. \"Doris: The Delft Object-oriented Radar Interferometric
00027  * software.\" In: proceedings 2nd ITC ORS symposium, August 1999. (cdrom)).
00028  *
00029  */
00030 /****************************************************************
00031  * $Source: /users/kampes/DEVELOP/DORIS/doris/src/RCS/matrixspecs.cc,v $
00032  * $Revision: 3.17 $
00033  * $Date: 2005/04/11 13:47:45 $
00034  * $Author: kampes $
00035  *
00036  * defenition of routines for matrices (helpers)
00037  ****************************************************************/
00038 
00039 
00040 #include "matrixbk.hh"
00041 #include <strstream>                    // 
00042 #include <algorithm>                    // max
00043 
00044 #ifdef __DEBUGMAT1
00045 #include <fstream>                      // type ofstream
00046 #endif
00047 
00048 //#if defined (__DEBUGMAT2) || defined (__DEBUGMAT1)
00049 #include <cstdlib>                      // exit
00050 //#endif
00051 
00052 
00053 
00054 // ______ message objects, global, set in main ______
00055 extern bk_messages matERROR;
00056 extern bk_messages matDEBUG;
00057 
00058 
00059 
00060 // ______ Prototype 1d fft in this file ______
00061 // ______ Only change this routine for other library for fft _____
00062 // ______ At the moment one can select VECLIB/FFTW/INTERNAL (slow?)
00063 // ______ VECLIB/FFTW call optimized routines for 2dfft's
00064 void four1(
00065         complr4 data[],
00066         int32 fftlength,
00067         int32 isign);
00068 
00069 // ______ Prototype veclib functions ______
00070 // ====== File matrixspecialization.c ======
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 //extern "C" { int32 s2dfft (real4*, real4*, int32*, int32*, int32*, int32*, int32*); }
00094 #endif
00095 
00096 // ______ USE FFTW IF AVAILABLE ______
00097 // first include complex then this to have complex type known?
00098 #ifdef __USE_FFTW_LIBRARY__
00099   #include <fftw3.h>                                  // fftw types and functions v3.0.1
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  * myispower2                                                   *
00118  *    Bert Kampes, 22-Mar-2000                                  *
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  *    matassert                                                 *
00132  * check if file is opened ok                                   *
00133  *    Bert Kampes, 16-Jun-1999                                  *
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   } // END matassert
00154 /****************************************************************
00155  *    matassert                                                 *
00156  * check if file is opened ok                                   *
00157  *    Bert Kampes, 16-Jun-1999                                  *
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   } // END matassert
00179 
00180 
00181 
00182 /****************************************************************
00183  *    fileci2tomatcr4                                           *
00184  * reads file format complex short (or 2*short)                 *
00185  * fills matrix complex real4                                   *
00186  *  window specified in system of file. win(10:20) means        *
00187  *  start at 10th line after 1st Byte on file.                  *
00188  * currentwindow should be passed to account for any offset     *
00189  * between what is linenumber of first line in file and byte = 1*
00190  *    Bert Kampes, 14-Jan-1999                                  *
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 // ______First account for possible offsets of file______
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   // ifstream ifile(file, ios::in | ios::ate | ios::binary);
00210   // g++ seems to haeve a prob. with this... BK 130700
00211   // ifstream ifile(file, ios::in | ios::app | ios::binary);
00212   //ifstream ifile(file, ios::in | ios::binary | ios::nocreate);
00213 #ifdef __NO_IOS_BINARY__
00214   ifstream ifile(file, ios::in);
00215 #else
00216   ifstream ifile(file, ios::in | ios::binary);
00217 #endif
00218 //  ifstream ifile;
00219 //  openfstream(ifile,file);
00220   ifile.seekg(0,ios::end);                      // pointer to end...
00221   matassert(ifile,file,__FILE__,__LINE__);
00222   const uint sizefile = ifile.tellg();          // opened ate
00223  
00224   if (sizefile==0)
00225     matERROR.print("filesize==0...");
00226 
00227 
00228 // ______Check input______
00229   const uint pixelsxbytes = sizefile/filelines;
00230   const uint SIZE = sizeof(compli16);           // size per element on disk (=4)
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   //const uint lines  = win.linehi - win.linelo + 1;
00236   //const uint pixels = win.pixhi - win.pixlo + 1;
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];                                          // pseudo complex type
00264 // ______Convert io_____
00265   for (int32 lin=0; lin<lines; ++lin)
00266     {
00267     ifile.seekg(start+filepixels*lin*sizepixel,ios::beg);    // data at row: lin
00268     for (int32 pix=0; pix<pixels; ++pix)
00269       {
00270       ifile.read((char*)el,SIZE);                       // element compli16 = 4
00271       Result(lin,pix)=complr4(el[0],el[1]);
00272       //Result(lin,pix)=complr4(fileci16.real(),fileci16.imag());
00273       }
00274     }
00275   ifile.close();
00276   } // END fileci2tomatcr4
00277 
00278 
00279 
00280 /****************************************************************
00281  *    magnitude                                                 *
00282  * Computes magnitude of complex matrix                         *
00283  *                                                              *
00284  *    Bert Kampes, 14-Jan-1999                                  *
00285  ****************************************************************/
00286 matrix<real4> magnitude(
00287         const matrix<complr4> &A)
00288   {
00289 #ifdef __DEBUGMAT2
00290   matDEBUG.print("magnitude.");
00291 #endif
00292 // _____Ensure stride one in memory______
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   } // END magnitude
00300 
00301 
00302 
00303 /****************************************************************
00304  *    intensity                                                 *
00305  * Computes intensity of complex matrix                         *
00306  *    Bert Kampes, 14-Jan-1999                                  *
00307  ****************************************************************/
00308 matrix<real4> intensity(
00309         const matrix<complr4> &A)
00310   {
00311 #ifdef __DEBUGMAT2
00312   matDEBUG.print("intensity.");
00313 #endif
00314 // _____Ensure stride one in memory______
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   } // END intensity
00322 
00323 
00324 
00325 #ifdef __USE_LAPACK_LIBRARY__
00326 /****************************************************************
00327  * choles(A);   cholesky factorisation lapack                   *
00328  * lower triangle of A is changed on output                     *
00329  * upper reamins un referenced                                  *
00330  *    Bert Kampes, 22-Feb-1999                                  *
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();                // leading dimension
00344   int32 n  = lda;                       // order of A
00345   int32 ierr;                           // 
00346   spotrf("U",&n,A[0],&lda,&ierr,1);     // Upper remains unchanged
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   } // END choles
00357 
00358 
00359 // ______ Else use simple intern algorithm ______
00360 // #elif defined (INTERNAL)
00361 #else
00362 /****************************************************************
00363  * choles(A);   cholesky factorisation internal implementation  *
00364  * lower triangle of A is changed on output                     *
00365  * upper reamins un referenced                                  *
00366  * this one is a lot slower then veclib and there may be more   *
00367  * efficient implementations.                                   *
00368  *    Bert Kampes, 11-Oct-1999                                  *
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   } // END choles internal, self bk
00398 #endif
00399 
00400 
00401 #ifdef __USE_LAPACK_LIBRARY__
00402 /****************************************************************
00403  * choles(A);   cholesky factorisation lapack                   *
00404  * lower triangle of A is changed on output                     *
00405  * upper reamins un referenced                                  *
00406  *    Bert Kampes, 22-Feb-1999                                  *
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();                // leading dimension
00420   int32 n  = lda;                       // order of A
00421   int32 ierr;                           // 
00422   dpotrf("U",&n,A[0],&lda,&ierr,1);     // Upper remains unchanged
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   } // END choles
00433 
00434 // ______ Else use simple intern algorithm ______
00435 // #elif defined (INTERNAL)
00436 #else
00437 /****************************************************************
00438  * choles(A);   cholesky factorisation internal implementation  *
00439  * lower triangle of A is changed on output                     *
00440  * upper reamins un referenced                                  *
00441  * this one is a lot slower then veclib and there may be more   *
00442  * efficient implementations.                                   *
00443  *    Bert Kampes, 11-Oct-1999                                  *
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   } // END choles internal, self bk
00473 #endif
00474 
00475 
00476 #ifdef __USE_LAPACK_LIBRARY__
00477 /****************************************************************
00478  * solvechol(A,rhs); solution of AX=rhs                         *
00479  *  cholesky factorisation lapack                               *
00480  * A contains cholesky factorisation of A                       *
00481  * rhs contains estimated X on output                           *
00482  *                                                              *
00483  * known BUG: niet mogelijk meer dan 1 rechterlid simultaan op  *
00484  *  te lossen, waarschijnlijk door fortran <> c tegenstelling   *
00485  *    Bert Kampes, 22-Feb-1999                                  *
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();                // leading dimension
00504   int32 ldb= B.lines();                 // leading dimension
00505   int32 nrhs= B.pixels();               // number of righthandsides
00506   int32 n  = lda;                       // order of A
00507   int32 ierr;                           // 
00508   spotrs("U",&n,&nrhs,A[0],&lda,B[0],&ldb,&ierr,1);     // factorisation stored in lower
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   } // END solvechol
00519 
00520 // ______ Else use simple (inaccurate, slow?) internal algorithm ______
00521 // #elif defined (INTERNAL)
00522 #else
00523 /****************************************************************
00524  * solvechol(A,rhs); solution of AX=rhs                         *
00525  *  cholesky factorisation internal implemetnation              *
00526  * A contains cholesky factorisation of A                       *
00527  * rhs contains estimated X on output                           *
00528  * there may be more efficient implementations.                 *
00529  *    Bert Kampes, 11-Oct-1999                                  *
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 // ______ Solve Ly=b, use B to store y ______
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 // ______ Solve Ux=y, use B to store unknowns ______
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   } // END solvechol
00570 #endif
00571 
00572 
00573 #ifdef __USE_LAPACK_LIBRARY__
00574 /****************************************************************
00575  * solvechol(A,rhs); solution of AX=rhs                         *
00576  *  cholesky factorisation lapack                               *
00577  * A contains cholesky factorisation of A                       *
00578  * rhs contains estimated X on output                           *
00579  *                                                              *
00580  * known BUG: niet mogelijk meer dan 1 rechterlid simultaan op  *
00581  *  te lossen, waarschijnlijk door fortran <> c tegenstelling   *
00582  *    Bert Kampes, 22-Feb-1999                                  *
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();              // leading dimension
00601   int32 ldb  = B.lines();               // leading dimension
00602   int32 nrhs = B.pixels();              // number of righthandsides
00603   int32 n    = lda;                     // order of A
00604   int32 ierr;                           // 
00605   dpotrs("U",&n,&nrhs,A[0],&lda,B[0],&ldb,&ierr,1);     // factorisation stored in lower
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   } // END solvechol
00616 
00617 // ______ Else use simple (inaccurate, slow?) internal algorithm ______
00618 // #elif defined (INTERNAL)
00619 #else
00620 /****************************************************************
00621  * solvechol(A,rhs); solution of AX=rhs                         *
00622  *  cholesky factorisation internal implemetnation              *
00623  * A contains cholesky factorisation of A                       *
00624  * rhs contains estimated X on output                           *
00625  * there may be more efficient implementations.                 *
00626  *    Bert Kampes, 11-Oct-1999                                  *
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 // ______ Solve Ly=b, use B to store y ______
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 // ______ Solve Ux=y, use B to store unknowns ______
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   } // END solvechol
00667 #endif
00668 
00669 
00670 #ifdef __USE_LAPACK_LIBRARY__
00671 /****************************************************************
00672  * invertchol(A); inversion by cholesky factorisation lapack    *
00673  * lower triangle of A is changed on output                     *
00674  * upper reamins un referenced                                  *
00675  *    Bert Kampes, 22-Feb-1999                                  *
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();                // leading dimension
00689   int32 n  = lda;                       // order of A
00690   int32 ierr;                           // 
00691   spotri("U",&n,A[0],&lda,&ierr,1);     // Upper remains unchanged
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   } // END invertchol
00702 
00703 
00704 // ______ Else use simple (inaccurate? slow?) internal algorithm ______
00705 // #elif defined (INTERNAL)
00706 #else
00707 /****************************************************************
00708  * invertchol(A); inversion by cholesky factorisation internal  *
00709  * lower triangle of A is changed on output                     *
00710  * upper remains un referenced                                  *
00711  * implementation with double loop probably can be improved easy*
00712  *    Bert Kampes, 11-Oct-1999                                  *
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 // ______ Compute inv(L) store in lower of A ______
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 // ______ Compute inv(A)=inv(LtL) store in lower of A ______
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];                 // transpose
00743         }
00744       A[j][i] = sum;
00745       }
00746     }
00747   } // END invertchol BK
00748 #endif
00749 
00750 
00751 #ifdef __USE_LAPACK_LIBRARY__
00752 /****************************************************************
00753  * invertchol(A); inversion by cholesky factorisation veclib    *
00754  * lower triangle of A is changed on output                     *
00755  * upper reamins un referenced                                  *
00756  *    Bert Kampes, 22-Feb-1999                                  *
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();                // leading dimension
00770   int32 n  = lda;                       // order of A
00771   int32 ierr;                           // 
00772   dpotri("U",&n,A[0],&lda,&ierr,1);     // Upper remains unchanged
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   } // END invertchol
00783 
00784 // ______ Else use simple (inaccurate? slow?) internal algorithm ______
00785 // #elif defined (INTERNAL)
00786 #else
00787 /****************************************************************
00788  * invertchol(A); inversion by cholesky factorisation internal  *
00789  * lower triangle of A is changed on output                     *
00790  * upper remains un referenced                                  *
00791  * implementation with double loop probably can be improved easy*
00792  *    Bert Kampes, 11-Oct-1999                                  *
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 // ______ Compute inv(L) store in lower of A ______
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 // ______ Compute inv(A)=inv(LtL) store in lower of A ______
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];                 // transpose
00823         }
00824       A[j][i] = sum;
00825       }
00826     }
00827   } // END invertchol BK
00828 #endif
00829 
00830 
00831 
00832 /****************************************************************
00833  *  A=norm(B)                                                   *
00834  *    Bert Kampes, 08-Feb-1999                                  *
00835  ****************************************************************/
00836 matrix<complr4> norm(
00837         const matrix<complr4> &A)
00838   {
00839 #ifdef __DEBUGMAT2
00840   matDEBUG.print("norm");
00841 #endif
00842 // ______Ensure stride one memory______
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   } // END norm
00850 
00851 
00852 
00853 /****************************************************************
00854  *    norm2                                                      *
00855  *    Bert Kampes, 08-Feb-1999                                  *
00856  ****************************************************************/
00857 real4 norm2(
00858         const matrix<complr4> &A)
00859   {
00860 #ifdef __DEBUGMAT2
00861   matDEBUG.print("norm2");
00862 #endif
00863 // ______Ensure stride one memory______
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   } // END norm2
00870 
00871 
00872 
00873 /****************************************************************
00874  *    norm2                                                      *
00875  *    Bert Kampes, 08-Feb-1999                                  *
00876  ****************************************************************/
00877 real4 norm2(
00878         const matrix<real4> &A)
00879   {
00880 #ifdef __DEBUGMAT2
00881   matDEBUG.print("norm2");
00882 #endif
00883 // ______Ensure stride one memory______
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   } // END norm2
00890 
00891 
00892 
00893 /****************************************************************
00894  * A = abs(B)                                                   *
00895  * Bert Kampes, 16-Feb-1999                                     *
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 // ______Ensure stride one memory______
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   } // END abs
00911 
00912 
00913 
00914 /****************************************************************
00915  * A = abs(B)                                                   *
00916  * Bert Kampes, 02-Jun-1999                                     *
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 // ______Ensure stride one memory______
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   } // END abs
00932 
00933 
00934 
00935 /****************************************************************
00936  * A=real(B)                                                    *
00937  *    Bert Kampes, 02-Mar-1999                                  *
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 // ______Ensure stride one memory______
00949   for (register int32 i=0; i<A.size(); i++)
00950     (*pntR++) = real(*pntA++);
00951   return Result;
00952   } // END real
00953 
00954 
00955 
00956 /****************************************************************
00957  * A=imag(B)                                                    *
00958  *    Bert Kampes, 02-Mar-1999                                  *
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 // ______Ensure stride one memory______
00970   for (register int32 i=0; i<A.size(); i++)
00971     (*pntR++) = imag(*pntA++);
00972   return Result;
00973   } // END imag
00974 
00975 
00976 
00977 /****************************************************************
00978  * A=angle(B)                                                   *
00979  *    Bert Kampes, 13-Apr-1999                                  *
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 // posible? imag(log(A)), matlab.
00991 // ______Ensure stride one memory______
00992   for (register int32 i=0; i<A.size(); i++)
00993     (*pntR++) = arg(*pntA++);
00994     //(*pntR++) = atan2((*pntA).imag(),(*pntA++).real());
00995   return Result;
00996   } // END angle
00997 
00998 
00999 
01000 /****************************************************************
01001  * A=angle2cmplx(B)                                             *
01002  *    Bert Kampes, 06-Oct-1999                                  *
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 // ______Ensure stride one memory______
01014   for (register int32 i=0; i<A.size(); i++)
01015     (*pntR++) = complr4(cos(*pntA),sin(*pntA++));
01016   return Result;
01017   } // END angle2cmplx
01018 
01019 
01020 
01021 /****************************************************************
01022  * dotmultconjphase(B,A)                                        *
01023  * subtract phase in real4 matrix A of complexr4 matrix B       *
01024  * by reference. (B .* conj(complr4(cos(phi,sin(phi))))         *
01025  *    Bert Kampes, 06-Oct-1999                                  *
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 // ______ Ensure stride one memory ______
01044 // --- looks a bit fishy: is pntA first increased from right to left,
01045 // --- or as i thought after sin is taken?
01046 // --- better would be to increase pntA in the for, together with i++
01047 // --- if we ever see a problem, it may be due to this?
01048   for (register int32 i=0; i<cint.size(); i++)
01049     (*pntC++) *= complr4(cos(*pntA),-sin(*pntA++));
01050   } // END dotmultconjphase
01051 
01052 
01053 
01054 /****************************************************************
01055  * A=coherence(complexinterferogram,norms,estwindowsizel,p)     *
01056  * returned is matrix size=(orig-winsizel,orig.pixels)          *
01057  * This routine is time consuming so no subroutines are called. *
01058  * a fft version is tried but does not speed up                 *
01059  *    Bert Kampes, 19-Apr-1999                                  *
01060  ****************************************************************/
01061 matrix<complr4> coherence(
01062         const matrix<complr4> &CINT, 
01063         const matrix<complr4> &NORMS,
01064         uint winL,
01065         uint winP)
01066   {
01067 // #define PLACE
01068 // #define FFTS
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   // ___ Allocate output matrix ___
01085   matrix<complr4> Result(CINT.lines()-winL+1,CINT.pixels());
01086   int32 leadingzeros  = (winP-1)/2;             // number of pixels=0 floor...
01087   int32 trailingzeros = (winP)/2;               // floor...
01088 
01089 // --- METHOD 1: --------------------------------------------- //
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 // --- end METHOD 1: ----------------------------------------- //
01114 
01115 
01116 
01117 // --- METHOD 2: --------------------------------------------- //
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     // ______ Compute sum over first block ______
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     // ______ Compute sum over rest of blocks ______
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 // --- end METHOD 2: ----------------------------------------- //
01153 
01154 
01155 
01156 // --- METHOD 3: --------------------------------------------- //
01157 // ====== Compute coherence by fft's ======
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;                       // power2, must be larger than winL
01164   int32 SIZEBLOCKP = 1024;                      // power2, must be larger than winP
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);                  // you should use theoretical expression.
01192 
01193   window win  = {0, sizeL-1, 0, sizeP-1};       // part from CINT,NORM
01194   window win2 = {0, sizeL-1, 0, sizeP-1};       // place where to put win
01195   register int32 i2 = 0;
01196   register int32 j2 = leadingzeros;
01197 
01198   for (;;) // ever
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 // ______ Convolution ______
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 // ______ Coherence ______
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 // ______ Update window ______
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;                                        // break loop
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;                            // this is last block
01296           win.linehi = CINT.lines() - 1;                // no resizing for now
01297           sizeL      = win.linehi - win.linelo + 1;
01298           }
01299         }
01300 
01301       else      // last block of line
01302         {
01303 cout << "BERT: newline==false\n";
01304         newline   = true;
01305         win.pixhi = CINT.pixels() - 1;                  // no resizing for now
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 // --- end METHOD 2: ----------------------------------------- //
01316 
01317 
01318   //cout << "Mean coherence: " << mean(Result); 
01319   return Result;
01320   } // END coherence
01321 
01322 
01323 
01324 /****************************************************************
01325  * A=coherence(complexinterferogram,norms,estwindowsizel,p)     *
01326  * returned is matrix size=(orig-winsizel,orig.pixels)          *
01327  * real4 coherence                                              *
01328  * This routine is time consuming so no subroutines are called. *
01329  *    Bert Kampes, 19-Apr-1999                                  *
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;             // number of pixels=0 floor...
01350   int32 trailingzeros = (winP)/2;               // floor...
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     // ______ Compute sum over first block ______
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     //Result(0,j) = norm(sum) / (power.real() * power.imag());
01370     real4 p = power.real() * power.imag();
01371     Result(0,j) = (p>0.0) ? sqrt(norm(sum)/p) : 0.0;
01372 
01373     // ______ Compute (relatively) sum over rest of blocks ______
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   } // END coherence2
01387 
01388 
01389 
01390 
01391 /****************************************************************
01392  * mysort2(A)                                                   *
01393  * sorts matrix first on column1, then on col2                  *
01394  * used to sort l,p,value matrices.                             *
01395  * calls to std c lib qsort, should better be c++ sort function?*
01396  * sorts in ascending order.                                    *
01397  * Bert Kampes, 10-Jan-2000                                     *
01398  ****************************************************************/
01399 int32 mycomp2 (const void *x1, const void *x2)
01400   {
01401   //  float *pntf = (float*) x1;
01402   //  cout << " BB : " << *pntf << " " << *(pntf+1) << endl;
01403   //  cout << " BB : " << *(float*)x1 << " " << *((float*)x1+1) << endl;
01404   int ret = 0;
01405   if      (*(float*)x1 < *(float*)x2 ) ret = -1; // first field smaller
01406   else if (*(float*)x1 > *(float*)x2 ) ret =  1;
01407   else if (*((float*)x1+1) < *((float*)x2+1) ) ret = -1; // x1=x2
01408   else if (*((float*)x1+1) > *((float*)x2+1) ) ret =  1; // x1=x2
01409   // else ret=0; // same x1,x2 and x1+1,x2+1
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   } // END mysort2
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   } // END mysort2
01436 
01437 
01438 
01439 //  /****************************************************************
01440 //   * mysort(A)                                                    *
01441 //   * sorts matrix on column col                                       *
01442 //   * calls to std c lib qsort, should better be c++ sort function?*
01443 //   * sorts in ascending order.                                    *
01444 //   * Bert Kampes, 10-Jan-2000                                     *
01445 //   ****************************************************************/
01446 //  int32 mycomp1 (const void *x1, const void *x2)
01447 //    {
01448 //  //  float *pntf = (float*) x1;
01449 //  //  cout << " BB : " << *pntf << " " << *(pntf+1) << endl;
01450 //    cout << " BB : " << *(float*)x1 << " " << *((float*)x1+1) << endl;
01451 //   
01452 //    int ret = 0;
01453 //    if      (*(float*)x1 < *(float*)x2 ) ret = -1; // first field smaller
01454 //    else if (*(float*)x1 > *(float*)x2 ) ret =  1;
01455 //    else if (*((float*)x1+columnnumber...) < *((float*)x2+columnnumber...) ) ret = -1; // x1=x2
01456 //    else if (*((float*)x1+columnnumber...) > *((float*)x2+columnnumber...) ) ret =  1; // x1=x2
01457 //    // else ret=0; // same x1,x2 and x1+1,x2+1
01458 //    return ret;
01459 //    }
01460 //  void mysort1 (matrix<real4> &A, column)
01461 //    {
01462 //  #ifdef __DEBUGMAT2
01463 //    matDEBUG.print("mysort2.");
01464 //  #endif
01465 //  #ifdef __DEBUGMAT1
01466 //    if (A.pixels() < 2)
01467 //      matERROR.print("mysort sorts only min. 2 cols.");
01468 //  #endif
01469 //    qsort(&A[0][0],A.lines(),A.pixels()*sizeof(real4),mycompare);
01470 //    } // END mysort1
01471 
01472 
01473 
01474 
01475 
01476 /****************************************************************
01477  * four1(complr4 *, length, isign)                              *
01478  * four1(&A[0][0], 128, 1)                                      *
01479  * either based on numerical recipes or veclib                  *
01480  * helper function for other fft routines, if no veclib         *
01481  * cooley-turkey, power 2, replaces input,                      *
01482  * isign=1: fft , isign=-1: ifft                                *
01483  * handling vectors should be simpler (lying, standing)         *
01484  * note that this is not a good implementation, only to get     *
01485  * doris software working without veclib.                       *
01486  *                                                              *
01487  * define SAMEASVECIB if you want the order of the coefficients *
01488  * of the fft the same as veclib. it seems this is not required *
01489  * for a good version of Doris, but in case of problems this    *
01490  * may be the solution.                                         *
01491  *                                                              *
01492  * VECLIB defines the FT same as matlab:                        * 
01493  *          N-1                                                 *
01494  *   X(k) = sum  x(n)*exp(-j*2*pi*k*n/N), 0 <= k <= N-1.        *
01495  *          n=0                                                 *
01496  *                                                              *
01497  * FFTW defines the same as Matlab, but inv. not normalized.    *
01498  * I don't know if the matrix must be allocated somehow, so for *
01499  * now we try only 1d ffts to build 2d too.                     *
01500  #%// BK 17-Sep-2003                                            *
01501  *                                                              *
01502  *    Bert Kampes, 12-Oct-1999                                  *
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 // --- common part to four1 for veclib,fftw, and internal ---
01510 void four1(
01511         complr4 data[],
01512         int32 fftlength,
01513         int32 isign)
01514   {
01515   // ______ Repair original routine ______
01516   //#ifdef __DEBUGMAT2 // called very often...
01517     //matDEBUG.print("four1: fft internal");
01518   //#endif
01519   #ifdef __DEBUGMAT1
01520   if (!myispower2(fftlength))
01521     matERROR.print("four1: length fft must be power of 2");// not really...
01522   if (abs(isign) != 1)
01523     matERROR.print("four1: isign should be (-)1.");
01524   #endif
01525 
01526 // ====== FFTW LIBRARY ROUTINE ======
01527 #ifdef __USE_FFTW_LIBRARY__
01528   // if no plan yet, create it depending on direction.
01529   // use guru interface for created plan on different data.
01530   // first check if dimension of old plan is still valid, or that the size of data changed,
01531   // and we thus have to create a new plan?
01532   // FFTW_ESTIMATE does not overwrite data during planning.
01533   // the plan can be stored in matrixclass as I do now, or as a static here...
01534   // that would mean only one plan, and no changes to matrixbk.h, but...
01535 // #define __NO_GURU_YET__
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);// deallocs, cannot use destroy when static?
01555 
01556 
01557   // === use guru interface to re-use a static plan ===
01558 #else
01559   if (isign == 1)
01560     {
01561     static int32 last_length_fwd = 0;
01562     static fftwf_plan fftwf_fwd_plan;// for fwd fft
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);// nerd-readable plan output
01579       matDEBUG.print(" ");
01580       #endif
01581       }
01582     // ___ Execute plan (guru interface) ___
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   // ___ Inverse transform ___
01593   else
01594     {
01595     static int32 last_length_inv = 0;
01596     static fftwf_plan fftwf_inv_plan;// for inverse fft
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);// nerd-readable plan output
01613       matDEBUG.print(" ");
01614       #endif
01615       }
01616     // ___ Execute plan (guru interface) ___
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     // ___ Inverse transform, normalization is not performed by fftw ___
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   }// END FOUR1
01633 #endif // guru interface for an unknown number of transforms
01634 #endif //fftw
01635 
01636 
01637 
01638   /**************************************************************
01639    * four1(complr4*, fftlength, iopt)                           *
01640    *    use veclib c1dfft                                       *
01641    *    data is complex vector (1,n)                            *
01642    *    isign=1: fft , isign=-1: ifft                           *
01643    **************************************************************/
01644 // ______ INTERNAL ROUTINE BASED ON NRC ______
01645 #ifndef __USE_FFTW_LIBRARY__
01646 #ifndef __USE_VECLIB_LIBRARY__
01647   register int32 i,j,n,mmax,m,istep;
01648   register complr4 bert;                // for SWAP function
01649   register double theta;                // trigionometric expansions
01650   register complr8 w,wp;                // trigionometric expansions
01651 
01652   #define SAMEASVECLIB
01653   #ifdef SAMEASVECLIB
01654   // bit slower, but for sure OK, spectral filtering etc. order may be important
01655   // to get same output as veclib, not sure if it makes a difference in computations,
01656   // but order of coefficients is defined differently
01657   // see also swap for isign=1.
01658   // order veclib: 0, -1/Nd .. 1/Nd, and this algorithm gives 0, 1/Nd .. -1/Nd
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   // ______ Check if input is standing or lying vector ______
01668   // index: j -> j/2 (floor) to correct for complex and first element
01669   // this should be repaired, made better later... within loop counter
01670   // define order differently
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   //#ifdef __GplusplusCOMPILER__                        // g++
01699   #ifdef __GNUC__
01700         bert       *= complr4(real4(w.real()),real4(w.imag())); // STUPID...
01701   #else                                         // assume aCC hp or so...
01702         bert       *= w; // not defined for g++ ???
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   // to get same output as veclib, not sure if it makes a difference in computations,
01714   // but order of coefficients is defined differently
01715   // this swap can be done better... above within loop somehow.
01716   // coefficients: 0 1 2 3 4 .. fftlength-1
01717   if (isign == 1)
01718     {for (i=1; i<fftlength/2; ++i)
01719       {SWAP(data[i],data[fftlength-i]);}}
01720   #endif
01721 
01722 
01723   // ______ Scale by fftlength if inverse transform ______
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   }// END FOUR1
01733 #endif //__USE_VECLIB_LIBRARY__
01734 #endif //__USE_FFTW_LIBRARY__
01735 
01736 
01737 
01738 // i prefer fftw for ffts
01739 #ifndef __USE_FFTW_LIBRARY__
01740 #ifdef __USE_VECLIB_LIBRARY__
01741   /**************************************************************
01742    * four1(complr4*, fftlength, iopt)                           *
01743    *    use veclib c1dfft                                       *
01744    *    data is complex vector (1,n)                            *
01745    *    isign=1: fft , isign=-1: ifft                           *
01746    *                                                            *
01747    *    If you want to include another library, intention is    *
01748    *    only to interchange this 1d fft routine.                *
01749    *                                                            *
01750    *    Bert Kampes, 22-Mar-2000                                *
01751    **************************************************************/
01752   int32 ierr = 0;
01753   static matrix<real4> work;                    // vector best lying
01754   if (int32(2.5*fftlength) != work.size())      // new size or first time??
01755     {
01756     work.resize(1,int32(2.5*fftlength));
01757     // ______ call initial work: iopt:=-3 ______
01758     int32 iopt = -3;                    // initialize work
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     // ______ Do actual 1d transform ____
01773   c1dfft(data,&fftlength,work[0],&isign,&ierr);
01774   if (ierr != 0) {cerr << "veclib: c1dfft: ierr = " << ierr << endl; exit(-1);}
01775   } // END four1
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 //#ifdef __USE_FFTW_LIBRARY__
01783 ///****************************************************************
01784 // * four1 with fftw plan given                                 *
01785 // #%// BK 17-Sep-2003                                          *
01786 // ****************************************************************/
01787 //void four1(
01788 //      complr4 data[],
01789 //        fftwf_plan fftwf_tmp_plan)
01790 //  {
01791 //  } // END four1
01792 //#endif
01793 
01794 
01795 
01796 
01797 /****************************************************************
01798  * fft(A,dim)                                                   *
01799  *    forward 1dfft over dim of A is returned in A by reference *
01800  *    if dim=1 fft is over all columns of A, if 2 over rows.    *
01801  *    data is stored major row order in memory, so dim=2 is     *
01802  *    probably much faster.                                     *
01803  *    fftlength should be power of 2                            *
01804  *    Bert Kampes, 22-Mar-2000                                  *
01805  ****************************************************************/
01806 void fft(matrix<complr4> &A, int32 dimension)
01807   {
01808   register int32 i;
01809   const int32 iopt = 1;                         // forward FFT
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         // note that fftw can be used directly to do this w/o data copying...
01821         matrix<complr4> VECTOR = A.getcolumn(i);// tmp copy, not very efficient maybe
01822         four1(&VECTOR[0][0],fftlength,iopt);// but generic.
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   } // END fft
01841 
01842 
01843 
01844 
01845 /****************************************************************
01846  * ifft(A,dim)                                                  *
01847  *    inverse 1dfft over dim of A is returned in A by reference *
01848  *    if dim=1 ifft is over all columns of A, if 2 over rows.   *
01849  *    data is stored major row order in memory, so dim=2 is     *
01850  *    probably much faster.                                     *
01851  *    fftlength should be power of 2                            *
01852  *    Bert Kampes, 22-Mar-2000                                  *
01853  ****************************************************************/
01854 void ifft(matrix<complr4> &A, int32 dimension)
01855   {
01856   register int32 i;
01857   const int32 iopt = -1;                        // inverse FFT (scaled)
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   } // END ifft
01888 
01889 
01890 
01891 // common part to fft2d for veclib, fftw, internal.
01892 /****************************************************************
01893  * fft2d(A);    2dfft veclib                                    *
01894  *    Bert Kampes, 01-Feb-1999                                  *
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;                        // leading dimension
01907   int32 iopt = 1;                       // forward transform if >= 0
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   } // END fft2d
01922 #endif
01923 #endif
01924 
01925 
01926 
01927 #ifdef __USE_FFTW_LIBRARY__ // prefer FFTW
01928   /**************************************************************
01929    * fft2d(A)                                                   *
01930    *    2dfft internal slow implementation,                     *
01931    *    four1 should be improved                                *
01932    *    optionally now fftw for 2d directly #%// BK 17-Sep-2003 *
01933    *    Bert Kampes, 12-Oct-1999                                *
01934    **************************************************************/
01935   // ______ fftw implementation 2d fft ______
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;// plan created each time, stupid.  test.
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);// nerd-readable plan output
01953   matDEBUG.print(" ");
01954   #endif
01955   // ___ execute plan ___
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);// also associated data?
01964   } // END fft2d
01965 #endif // select FFTW
01966 
01967 
01968 // ====== Internal implementation using series of 1d fft-s ======
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);                             // forward transform over rows
01975   fft(A,1);                             // forward transform over columns
01976   } // END fft2d
01977 #endif // select internal
01978 #endif // select internal
01979 
01980 
01981 
01982 
01983 // common part to ifft2d for veclib, fftw, internal.
01984 /****************************************************************
01985  * ifft2d(A);   2dfft veclib                                    *
01986  *    Bert Kampes, 01-Feb-1999                                  *
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;                        // leading dimension
01999   int32 iopt = -1;                      // invers transform if < 0
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   } // END ifft2d
02014 #endif
02015 #endif
02016 
02017 
02018 
02019 #ifdef __USE_FFTW_LIBRARY__   // prefer fftw
02020 // ====== use internal slow implementation ======
02021 /****************************************************************
02022  * ifft2d(A);   2dfft internal slow implementation,             *
02023  * should be improved                                           *
02024  *    Bert Kampes, 12-Oct-1999                                  *
02025  ****************************************************************/
02026   // ______ fftw implementation 2d fft ______
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;// plan created each time, stupid.  test.
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);// nerd-readable plan output
02044   matDEBUG.print(" ");
02045   #endif
02046   // ___ execute plan ___
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);// also associated data?
02055   // ___ normalize inverse transform, not done by fftw ___
02056   #ifdef __DEBUGMAT2
02057   matDEBUG.print("BK: FFTW: normalizing inverse transform."); 
02058   #endif
02059   A /= complr4(nx*ny);
02060   //for (int32 i=0; i<ny; ++i)
02061   //  for (int32 j=0; i<nx; ++j)
02062   //#ifdef __GNUC__
02063   //        A[i][j] /= nx*ny;
02064   //#else
02065   //        A[i][j] /= complr4(nx*ny);
02066   //#endif
02067   } // END ifft2d
02068 #endif
02069 
02070 
02071 
02072 // ====== Internal implementation using series of 1d fft-s ======
02073 #ifndef __USE_FFTW_LIBRARY__
02074 #ifndef __USE_VECLIB_LIBRARY__
02075   #ifdef __DEBUGMAT2
02076     matDEBUG.print("ifft2d: internal");
02077   #endif
02078   ifft(A,2);                            // inverse transform over rows
02079   ifft(A,1);                            // inverse transform over columns
02080   } // END ifft2d
02081 #endif
02082 #endif
02083 
02084 
02085 
02086 /****************************************************************
02087  * B=oversample(A, factor);             harmonic interpolation  *
02088  *    Bert Kampes, 01-Feb-1999                                  *
02089  *    removed from code.                                        *
02090  * B=oversample(A, factorrow, factorcol);                       *
02091  *    2 factors possible, extrapolation at end.                 *
02092  *    no vectors possible.                                      *
02093  *    Bert Kampes, 28-Mar-2000                                  *
02094  ****************************************************************/
02095 matrix<complr4> oversample(
02096         matrix<complr4> A,              // not by reference, changed by fft
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;       // numrows of output matrix
02108   const uint P2    = factorcol*p;       // columns of output matrix
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);                           // 1d fourier transform per row
02125     for (i=0; i<l; ++i)                 // divide by 2 'cause even fftlength
02126 #ifdef __GNUC__
02127       A(i,halfp) *= half;               // complex *=?
02128 #else
02129       A(i,halfp) *= complr4(half);              // not tested...complex *=?
02130 #endif
02131     const window winA1(0, l-1, 0, halfp);       // zero padding windows
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);                        // inverse fft per row
02137     }
02138   else if (factorcol==1)
02139     {
02140     fft(A,1);                           // 1d fourier transform per column
02141     for (i=0; i<p; ++i)                 // divide by 2 'cause even fftlength
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);       // zero padding windows
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);                        // inverse fft per row
02153     }
02154   else
02155     {
02156     fft2d(A);                           // A=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);   // zero padding windows
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   //return (factorrow*factorcol)*Res;           // overhead?
02186   } // END oversample
02187 
02188 
02189 
02190 /****************************************************************
02191  * B=oversample(A, f1, f2);     calls complex oversampling      *
02192  *    Bert Kampes, 28-Mar-2000                                  *
02193  ****************************************************************/
02194 matrix<real4> oversample(
02195         const matrix<real4> &A,
02196         uint factorrow,
02197         uint factorcol)
02198   {
02199   // not efficient, mat2cr4 allocates new, and oversample copies this.
02200   matrix<complr4> TMP = oversample(mat2cr4(A),factorrow,factorcol);
02201   return real(TMP);                     // imag == 0
02202   }
02203 
02204 
02205 
02206 /****************************************************************
02207  * dotmult(complr4* pnt, const matrix<real4>&B, stride)         *
02208  * multiply memory pointed to by content of B, with mem. stride *
02209  *    Bert Kampes, 04-Apr-2000                                  *
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 //#ifdef __DEBUGMAT1
02220 //  if (pntR+B.size()*stride>??)
02221 //    matERROR.print("to large....");
02222 //#endif
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   } // dotmult
02234 
02235 
02236 
02237 // ====== First declare (define?) specialized function ======
02238 // ______ Before template function ______
02239 #ifdef __USE_VECLIB_LIBRARY__
02240 /****************************************************************
02241  * C = A * B; real4                                             *
02242  * Bert Kampes, 14-Jan-1999                                     *
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   } // END *
02264 #endif
02265 
02266 
02267 #ifdef __USE_VECLIB_LIBRARY__
02268 /****************************************************************
02269  * C = A * B; real8                                             *
02270  * Bert Kampes, 14-Jan-1999                                     *
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   } // END *
02292 #endif
02293 
02294 
02295 #ifdef __USE_VECLIB_LIBRARY__
02296 /****************************************************************
02297  * C = A * B; complex real4                                             *
02298  * Bert Kampes, 14-Jan-1999                                     *
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   } // END *
02320 #endif
02321 
02322 
02323 #ifdef __USE_VECLIB_LIBRARY__
02324 /****************************************************************
02325  * C = A * B; complex real8                                             *
02326  * Bert Kampes, 14-Jan-1999                                     *
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   } // END *
02348 #endif
02349 
02350 
02351 
02352 #ifdef __USE_VECLIB_LIBRARY__
02353 /****************************************************************
02354  * C=matTxmat(A,B) C=trans(A)*B; specialized for veclib         *
02355  *    Bert Kampes, 22-Feb-1999                                  *
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   } // END matTxmat
02375 #endif
02376 
02377 
02378 
02379 #ifdef __USE_VECLIB_LIBRARY__
02380 /****************************************************************
02381  * C=matTxmat(A,B) C=trans(A)*B; specialized for veclib         *
02382  *    Bert Kampes, 22-Feb-1999                                  *
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   } // END matTxmat
02402 #endif
02403 
02404 
02405 
02406 #ifdef __USE_VECLIB_LIBRARY__
02407 /****************************************************************
02408  * C=matTxmat(A,B) C=trans(A)*B; specialized for veclib         *
02409  *    Bert Kampes, 22-Feb-1999                                  *
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   } // END matTxmat
02429 #endif
02430 
02431 
02432 
02433 #ifdef __USE_VECLIB_LIBRARY__
02434 /****************************************************************
02435  * C=matTxmat(A,B) C=trans(A)*B; specialized for veclib         *
02436  *    Bert Kampes, 22-Feb-1999                                  *
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   } // END matTxmat
02456 #endif
02457 
02458 
02459 
02460 #ifdef __USE_VECLIB_LIBRARY__
02461 /****************************************************************
02462  * C=matxmatT(A,B) C=A*trans(B); specialized for veclib         *
02463  *    Bert Kampes, 22-Feb-1999                                  *
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   } // END matxmatT
02483 #endif
02484 
02485 
02486 
02487 #ifdef __USE_VECLIB_LIBRARY__
02488 /****************************************************************
02489  * C=matxmatT(A,B) C=A*trans(B); specialized for veclib         *
02490  *    Bert Kampes, 22-Feb-1999                                  *
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   } // END matxmatT
02510 #endif
02511 
02512 
02513 
02514 #ifdef __USE_VECLIB_LIBRARY__
02515 /****************************************************************
02516  * C=matxmatT(A,B) C=A*trans(B); specialized for veclib         *
02517  *    Bert Kampes, 22-Feb-1999                                  *
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   } // END matxmatT
02537 #endif
02538 
02539 
02540 
02541 #ifdef __USE_VECLIB_LIBRARY__
02542 /****************************************************************
02543  * C=matxmatT(A,B) C=A*trans(B); specialized for veclib         *
02544  *    Bert Kampes, 22-Feb-1999                                  *
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   } // END matxmatT
02564 #endif
02565 
02566 
02567 
02568 /****************************************************************
02569  * C *= A;      pointwise multiplication, CR4 R4                *
02570  #%// BK 26-Oct-2000
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   } // END *= matrices CR4 with R4
02587 
02588 
02589 
02590 /****************************************************************
02591  * C /= A;      pointwise multiplication, CR4 R4                *
02592  #%// BK 27-Oct-2000                                            *
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   } // END /= matrices CR4 with R4
02609 
02610 
02611 
02612 /****************************************************************
02613  * C += A;      pointwise multiplication, CR4 R4                *
02614  #%// BK 27-Oct-2000                                            *
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   } // END += matrices CR4 with R4
02631 
02632 
02633 
02634 /****************************************************************
02635  * C -= A;      pointwise subtraction, CR4 R4           *
02636  #%// BK 27-Oct-2000                                            *
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   } // END -= matrices CR4 with R4
02653 
02654 
02655 
02656 /****************************************************************
02657  * C=diagxmat(vec,B) C=diag(vec) * B; (diag R4, mat CR4)        *
02658  *    Bert Kampes, 22-Feb-1999                                  *
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)        // standing
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         //Result.data[i][j] *= diag.data[i][0];
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         //Result.data[i][j] *= diag.data[0][i];
02686     }
02687   return Result;
02688   } // END diagxmat
02689 
02690 
02691 
02692 /****************************************************************
02693  * A=cos(B)                                                     *
02694  #%// BK 09-Nov-2000                                            *
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   // ______Ensure stride one memory______
02706   for (register int32 i=0; i<A.size(); i++)
02707     (*pntR++) = cos(*pntA++);
02708   return Result;
02709   } // END cos
02710 
02711 
02712 
02713 /****************************************************************
02714  * A=cos(B)                                                     *
02715  #%// BK 09-Nov-2000                                            *
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   // ______Ensure stride one memory______
02727   for (register int32 i=0; i<A.size(); i++)
02728     (*pntR++) = cos(*pntA++);
02729   return Result;
02730   } // END cos
02731 
02732 
02733 
02734 /****************************************************************
02735  * A=sin(B)                                                     *
02736  #%// BK 09-Nov-2000                                            *
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   // ______Ensure stride one memory______
02748   for (register int32 i=0; i<A.size(); i++)
02749     (*pntR++) = sin(*pntA++);
02750   return Result;
02751   } // END sin
02752 
02753 
02754 
02755 /****************************************************************
02756  * A=sin(B)                                                     *
02757  #%// BK 09-Nov-2000                                            *
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   // ______Ensure stride one memory______
02769   for (register int32 i=0; i<A.size(); i++)
02770     (*pntR++) = sin(*pntA++);
02771   return Result;
02772   } // END sin
02773 
02774 
02775 
02776 /****************************************************************
02777  * B=mat2cr4(A); conversion to complex<real4>                   *
02778  *    Bert Kampes, 01-Feb-1999                                  *
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   } // END mat2cr4 (real4)
02793 
02794 
02795 
02796 /****************************************************************
02797  * A=mat2cr4(B,C)                                               *
02798  #%// BK 09-Nov-2000                                            *
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   } // END mat2cr4 complr4 (real4,real4)
02819 
02820 
02821 
02822 /****************************************************************
02823  * A=mat2cr4(r8B,r8C)                                           *
02824  #%// BK 09-Nov-2000                                            *
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   } // END mat2cr4(real8,real8)
02845 
02846 
02847 
02848 //  /****************************************************************
02849 //   * A=mat2cr4(B,C)                                                   *
02850 //   #%// BK 09-Nov-2000                                                *
02851 //   ****************************************************************/
02852 //  matrix<complr8> mat2cr8(
02853 //      const matrix<real8>& A,
02854 //      const matrix<real8>& B)
02855 //    {
02856 //    #ifdef __DEBUGMAT1
02857 //    if (A.lines()!=B.lines() || A.lines()!=B.lines())
02858 //      matERROR.print("operator complr8, input matrices not same size");
02859 //    #endif
02860 //    #ifdef __DEBUGMAT2
02861 //    matDEBUG.print("operator complr8(r8,r8)");
02862 //    #endif
02863 //    matrix<complr8> Result(A.lines(),A.pixels());
02864 //    real8   *pntA = A[0];
02865 //    real8   *pntB = B[0];
02866 //    complr8 *pntR = Result[0];
02867 //    for (register int32 i=0; i<Result.size(); i++)
02868 //      (*pntR++) = complr8(*pntA++,*pntB++);
02869 //    return Result;
02870 //    } // END mat2cr8 (real8,r8)
02871 
02872 
02873 
02874 
02875 // ============================================================
02876 // explicit instantitiate members not required...
02877 //template class matrix<float>;
02878 //template class matrix<complr4>;
02879 
02880 // explicit instantitiate friends...
02881 //template <class Type>
02882 //  matrix<Type>  operator * (const matrix<Type>& A, const matrix<Type>& B);
02883 
02884 //template matrix<complr4>  operator * <complr4> (const matrix<complr4>& A, const matrix<complr4>& B);
02885 
02886 
02887 // ok
02888 //template matrix<complr4>  operator * (const matrix<complr4>& A, const matrix<complr4>& B);

Generated on Fri Apr 22 15:57:56 2005 for Doris by doxygen 1.3.6