EMAN2
utilnum.cpp
Go to the documentation of this file.
00001 #include <math.h>
00002 #include <stdlib.h>
00003 #include <stdio.h>
00004 #include <math.h>
00005 #include <string.h>
00006 #include <ctype.h>
00007 #include <float.h>
00008 
00009 #include "lapackblas.h"
00010 #include "utilnum.h"
00011 
00012 int asta2(float *img, int nx, int ny, int ri, double *abaloc, int *klploc)
00013 {
00014     // calculate the sum of the background pixels, and then set the intensity
00015     // of these pixels to zero. Also count the number of background pixels.
00016     // A background pixel is a pixel outside of the circle with radius ri
00017     // This is done for images assigned to a processor
00018     //
00019     int xcent = (nx / 2) + 1;
00020     int ycent = (ny / 2) + 1;
00021     int r_squared = ri*ri;
00022 
00023     int x_summand, y_summand;
00024 
00025     for ( int i = 0 ; i < nx ; ++i ) {
00026         x_summand = (i-xcent) * (i-xcent);
00027         for ( int j = 0 ; j < ny ; ++j ) {
00028             y_summand = (j-ycent) * (j-ycent);
00029             if ( x_summand + y_summand > r_squared ) {
00030                 *abaloc += (double) img[j*nx + i];
00031                 //chao set the background to zero
00032                 img[j*nx+i]=0.0;
00033                 ++*klploc;
00034             }
00035         }
00036     }
00037 
00038     return 0;
00039 }
00040 
00041 int ifix(float a)
00042 {
00043     int ia = 0;
00044     if (a >= 0.0) {
00045        ia = (int)floor(a);
00046     }
00047     else {
00048        ia = (int)ceil(a);
00049     }
00050     return ia;
00051 }
00052 
00053 // ================= imported from direct_methods.c =================
00054 
00055 /* direct_methods.c
00056  * ------------------------------------------------------------------
00057  *  This file contains SVD based methods to solve discrete
00058  *  ill-posed problems.
00059  *
00060  *------------------------------------------------------------------
00061  */
00062 
00063 
00064 /*-------SVD--------------------------------------------------------*/
00065 void svd(float **A, int m, int n, float **U, float *s, float **V)
00066 /*
00067  * Computes the SVD of an m-by-n matrix A:
00068  *
00069  *        A = USV^T
00070  *
00071  *  NOTE: We assume m >= n
00072  *
00073  *  Input:  A    : m-by-n matrix
00074  *          m, n : row and column dimensions of A
00075  *
00076  *  Output: U    : Orthogonal matrix containing left singular vectors.
00077  *          s    : vector containing singular values of A.
00078  *               : That is, s = diag(S)
00079  *          V    : Orthogonal matrix containing right singular vectors.
00080  */
00081  
00082 { 
00083   char jobu = 'A';
00084   char jobv = 'A';
00085   
00086   int i;
00087 
00088 // chao something funny about min/max here... hardwire lwork to 100 for now  
00089 //  integer lwork = max(3*min(m,n) + max(m,n), 5*min(m,n));
00090   integer lwork = 100;
00091   integer lda=m; 
00092   integer ldu=m; 
00093   integer ldv=n;
00094   integer info;
00095   integer mm = m;
00096   integer nn = n;  
00097 
00098   float ** A_column = matrix(n,m);
00099   float * work;
00100   
00101   work = myvector((int) lwork);
00102 
00103   /* We have to column orient the matrix going in because dgesvd expects the 
00104    * matrix to be column oriented and in C, the matrix is row oriented. */
00105  column_orient(A, m, n, A_column);
00106   
00107  /* The LAPACK routines destroy the input matrix, so use a
00108   * copy of it.
00109   *
00110   * Check to see if we are using single or float precision,
00111   * and use appropriate LAPACK routine to compute SVD.
00112   */
00113 
00114  sgesvd_(&jobu, &jobv, &mm, &nn, *A_column, &lda, s, *U, &ldu, *V, &ldv, 
00115          work,  &lwork, &info);
00116 
00117 //sgesvd_(jobu, jobv, mm, nn, *A_column, lda, s, *U, ldu, *V, ldv, &info);
00118 
00119   /* We have to transpose the U matrix to make it row oriented again. 
00120    * We don't worry about V because we would have had to transpose it anyways. */
00121         transpose_matrix(U, m, m);
00122   
00123   free_vector(work);
00124   free_matrix(A_column);
00125 
00126   
00127 }
00128 
00129 /*-------TSVD--------------------------------------------------------*/
00130 void tsvd(float **U, float *s, float **V, float *b, int m, int n, float omega, float *x, float *tol)
00131    /*
00132     *  Computes the TSVD solution
00133     *
00134     *  NOTE: We assume m >= n
00135     *
00136     *  Input:  U    : Orthogonal matrix containing left singular vectors.
00137     *          s    : vector containing singular values of A.
00138     *               : That is, s = diag(S)
00139     *          V    : Orthogonal matrix containing right singular vectors.
00140     *          b    : Right hand side vector
00141     *          m, n : row and column dimensions of A.
00142     *          tol  : Truncation tolerance.  If tol < 0, then
00143     *                   GCV is used to choose it.
00144     *
00145     *  Output: x    : TSVD solution
00146     *          tol  : Truncation tolerance.
00147     */
00148 {    
00149     int i, k;
00150     float * sinv = myvector(n);
00151     float * bhat = myvector(m);
00152 
00153     matvec_multT(U, b, m, m, bhat);
00154 
00155 //printf("m=%d, n=%d\n", m,n);
00156      /*
00157       * Use GCV to choose tol.
00158       */
00159     k = 0;
00160     if (*tol < (float) 0.0)
00161       gcv_tsvd(s, bhat, m, n, omega, tol);
00162 
00163     /*
00164      * Now we need to compute the the filtered inverse of the
00165      * singular values.  That is, 1/s(i) if s(i) >= tol, and
00166      * 0 otherwise.  I'm not sure this is the best way to do this
00167      * in Fortran -- it would be better to not have to check each
00168      * time through the loop.  But, I do not want to assume the
00169      * the singular values are ordered -- they may not be in 
00170      * some situations. Maybe a pre-sort might be helpful???
00171      */
00172  
00173     for (i=0; i<n; i++){
00174       if (fabs(s[i]) >= *tol){
00175         sinv[i] = (float) 1.0 / s[i];
00176         k++;
00177       }
00178     }
00179 *tol = k;
00180 //printf("truncindex = %f\n", *tol);
00181 //print_vector(s,n);
00182 
00183     /*
00184      * Now the TSVD solution is easy to compute.
00185      */
00186     for (i=0; i<n; i++)
00187       bhat[i] = sinv[i]*bhat[i];
00188     
00189     matvec_mult(V, bhat, n, n, x);
00190     
00191     
00192     free_vector(sinv);
00193     free_vector(bhat);
00194 }
00195     
00196     
00197     
00198 /*-----GCV_TSVD------------------------------------------------*/
00199 
00200 void gcv_tsvd(float *s, float *bhat, int m, int n, float omega, float * tol)
00201    /*
00202     *  This function uses the GCV function to choose a default truncation tolerance for TSVD regularization
00203     *
00204     *  Use GCV to compute default regularization parameter.
00205     *
00206     * Input:    s : Vector containing singular values.
00207     *        bhat : Vector containing spectral coefficients of right
00208     *                 hand side vector.  That is, bhat = U^T*b.
00209     *
00210     * Output: tol : Scalar truncation tolerance for TSVD.
00211     */
00212 {
00213     float * bhat_copy = myvector(m);
00214     float * rho = myvector(m-1);
00215     float * G = myvector(m-1);
00216     int i, kk=0;
00217     float ni;
00218     /*
00219      *  NEED TO: Make sure singular values are sorted from
00220      *           largest to smallest, and do the same for the
00221      *           bhat values.  In MATLAB this looks like:
00222      *              [s, idx] = sort(abs(s)); 
00223      *              s = flipud(s); 
00224      *              idx = flipud(idx);
00225      *              bhat = abs( bhat(idx) );
00226      *
00227      * Initialize residual and GCV function values.
00228      */
00229     
00230     copy_vector(bhat, bhat_copy, m);
00231     for (i=0; i<m; i++)
00232       bhat_copy[i] = fabs( bhat[i] );
00233     
00234     rho[m-2] = bhat_copy[m-1]*bhat_copy[m-1];
00235 
00236     /*
00237      * Recursively compute residual and GCV function values.
00238      */
00239     for (i = m-3; i>-1; i--){
00240       rho[i] = rho[i+1] + bhat_copy[i+1]*bhat_copy[i+1];
00241     }
00242 
00243 //    print_vector(bhat_copy, m-1);
00244     for(i=0; i< m-1; i++){
00245       ni = (float) m - omega*(i+1);
00246       G[i] = (float) n* rho[i]/(ni*ni);
00247     }
00248  // print_vector(G,m-1);
00249     /*
00250      * Now find minimum of the discrete GCV function.
00251      */
00252     for (i=0; i < m-2; i++){
00253       if (G[i+1]<G[i])
00254         kk = i+1;
00255     }
00256 
00257     *tol = s[kk];
00258  // printf("kk = %d, s[kk] = %f\n", kk, *tol);
00259 
00260     free_vector(bhat_copy);
00261     free_vector(rho);
00262     free_vector(G);
00263 }
00264 
00265 
00266 /*------TIKHONOV------------------------------------------------------------*/
00267 
00268 void tikhonov(float **U, float *s, float **V, float *b, int m, int n, float omega, float *x, float * alpha)
00269   /*
00270    *  Computes the Tikhonov solution
00271    *
00272    *  NOTE: We assume m >= n
00273    *
00274    *  Input:  U    : Orthogonal matrix containing left singular vectors.
00275    *          s    : vector containing singular values of A.
00276    *               : That is, s = diag(S)
00277    *          V    : Orthogonal matrix containing right singular vectors.
00278    *          b    : Right hand side vector
00279    *          m, n : row and column dimensions of A.
00280    *        alpha  : Regularization parameter.  If alpha < 0, then
00281    *                   GCV is used to choose it.
00282    *        omega  : Omega parameter
00283    *
00284    *  Output: x    : Tikhonov solution
00285    *        alpha  : Regularization parameter
00286    */
00287 {
00288     int i, MaxIts = 500;
00289     float * bhat = myvector(m);
00290     float * xhat = myvector(n);
00291     float * D = myvector(n);
00292 
00293     matvec_multT(U, b, m, m, bhat);
00294 
00295 /*printf("b:\n");
00296     print_vector(b, m);
00297     printf("m = %d, n=%d, omega=%f\n", m,n,omega);
00298 printf("bhat:\n");
00299     print_vector(bhat, m);
00300 printf("s:\n");
00301     print_vector(s, n);
00302 */
00303 
00304 
00305      /*
00306       * Use GCV to choose alpha.
00307       */
00308     
00309     if (*alpha < (float)0.0)
00310        *alpha = GCVmin_Tik(s, bhat, m, n, MaxIts, omega);
00311 
00312    /*
00313     * Now we need to compute the the filtered Tikhonov solution
00314     */
00315     for (i=0; i<n;i++){
00316      D[i] = fabs(s[i])*fabs(s[i]) + (*alpha)*(*alpha);
00317      xhat[i] = (s[i] * bhat[i])/D[i];
00318     }
00319  
00320     matvec_mult(V, xhat, n, n, x);
00321 
00322     free_vector(bhat);
00323     free_vector(xhat);
00324     free_vector(D);
00325 }
00326 
00327 /*------GCVmin_Tik------------------------------------------------------------*/
00328 
00329 float GCVmin_Tik(float *s, float *bhat, int m, int n, int MaxIts, float omega)
00330    /*
00331     *
00332     *  Use a combination of Golden Section Search and Parabolic
00333     *  Interpolation to compute the minimum of the GCV 
00334     *  function to find a default regularization parameter (see MATLAB's fminbnd)
00335     *
00336     * Input:    s : Vector containing singular or spectral values of A.
00337     *        bhat : Vector containing spectral coefficients of right
00338     *                 hand side vector.  That is, bhat = U^T*b.
00339     *       MaxIts : Maximum number of iterations to find the minimum
00340     *
00341     * Output: xf : minimum of the GCV function for Tikhonov.
00342     */
00343 {
00344     float ax, bx, seps, c, a, b, v, w, sf, d, e, x, fx, fv, fw, xm, tol, tol1, tol2, r, q, p, si, fu, fval, xf;
00345     int funccount, iter, maxfun, gs;
00346     
00347 /*float tt;
00348 
00349 tt = GCVfun_Tik(0.01, s, bhat, m,n, omega);
00350 printf("alpha = 0.01;  GCV(a) = %f\n", tt);
00351 tt = GCVfun_Tik(0.59, s, bhat, m,n, omega);
00352 printf("alpha = 0.59;  GCV(a) = %f\n", tt);
00353 tt = GCVfun_Tik(.999, s, bhat, m,n, omega);
00354 printf("alpha = .999;  GCV(a) = %f\n", tt);
00355 tt = GCVfun_Tik(1.2, s, bhat, m,n, omega);
00356 printf("alpha = 1.2;  GCV(a) = %f\n", tt);
00357 */
00358 
00359     ax = (float) 0.0;
00360  //   bx = (float)1.0;
00361    bx = s[0];
00362 
00363     /*  Initialize some values */
00364     seps = sqrt(FLT_EPSILON);
00365     c = (float)0.5*((float)3.0 - sqrt((float)5.0));
00366     a = ax;
00367     b = bx;
00368     v = a + c*(b-a);
00369     w = v;
00370     xf = v;
00371     d = (float)0.0;
00372     e = (float)0.0;
00373     x = xf; 
00374 
00375     fx = GCVfun_Tik(x, s, bhat, m,n, omega);
00376 
00377     fv = fx;
00378     fw = fx;
00379     xm = (float)0.5*(a+b);
00380     tol = (float).0001;
00381     tol1 = seps*fabs(xf) + tol/(float)3.0;
00382     tol2 = (float)2.0*tol1;
00383     funccount = 0;
00384     iter = 0;
00385     maxfun = 500;
00386 
00387     for (iter = 1; iter< MaxIts; iter++){
00388       gs = 1;
00389       /* Is a parabolic fit possible*/
00390       if (fabs(e) > tol1){
00391         gs = 0;
00392         r = (xf-w)*(fx-fv);
00393         q = (xf-v)*(fx-fw);
00394         p = (xf-v)*q-(xf-w)*r;
00395         q = (float)2.0*(q-r);
00396         if (q > (float)0.0)
00397           p = -p;
00398         q = fabs(q);
00399         r = e;
00400         e = d;
00401         /*Is the parabola acceptable*/
00402         if ( (fabs(p)<fabs((float)0.5*q*r)) && (p>q*(a-xf)) && (p<q*(b-xf)) ){
00403           /* Yes, parabolic interpolation step*/
00404           d = p/q;
00405           x = xf+d;
00406           if (((x-a) < tol2) || ((b-x) < tol2)){
00407             if(xm-xf <0)
00408               si = -(float)1.0;
00409             else
00410               si = (float)1.0;
00411             
00412             d = tol1*si;
00413           }
00414         }
00415         else{
00416           /* Not acceptable, must do a golden section step*/
00417           gs = 1;
00418         }
00419       }
00420       if (gs){
00421         /* A golden-section step is required*/
00422         if (xf >= xm)
00423           e = a-xf;
00424         else
00425           e = b-xf;
00426 
00427         d = c*e;
00428       }
00429       /* The function must not be evaluated too close to xf*/
00430      if(d <0)
00431        si = -(float)1.0;
00432      else
00433        si = (float)1.0;
00434      
00435       x = xf + si * max( fabs(d), tol1 );
00436       fu = GCVfun_Tik( x, s, bhat, m,n, omega);
00437       funccount = funccount + 1;
00438       if (fu <= fx){
00439         if (x >= xf)
00440           a = xf;
00441         else
00442           b = xf;
00443    
00444         v = w;
00445         fv = fw;
00446         w = xf;
00447         fw = fx;
00448         xf = x;
00449         fx = fu;
00450     }
00451       else{
00452         if (x < xf)
00453           a = x;
00454         else 
00455           b = x;
00456         
00457         if ( (fu <= fw) || (w == xf) ){
00458           v = w;
00459           fv = fw;
00460           w = x;
00461           fw = fu;
00462         }
00463         else if ((fu <= fv) || (v == xf) || (v == w)){
00464           v = x;
00465           fv = fu;
00466         }
00467       }
00468       xm = (float)0.5*(a+b);
00469       tol1 = seps*fabs(xf) + tol/(float)3.0;
00470       tol2 = (float)2.0*tol1;
00471       if (funccount >= maxfun || fabs(xf-xm) <= (tol2 - (float)0.5*(b-a))){
00472         fval = fx;
00473         return xf;
00474       }
00475 }
00476     fval = fx;
00477     return xf;
00478 }
00479 
00480 
00481 /*---------GCVFUN_TIK---------------------------------------------------------*/
00482 
00483 float GCVfun_Tik(float alpha, float * s, float * bhat, int m, int n, float omega)
00484     /*
00485      *   This function evaluates the GCV function for Tikhonov regularization.
00486      *   Assume the problem has the form b = Ax + noise, with [U, S, V] = svd(A).
00487      *
00488      * Input: alpha: value at which to evaluate the GCV function
00489      *           s : Vector containing singular or spectral values of A.
00490      *        bhat : Vector containing spectral coefficients of right
00491      *                 hand side vector.  That is, bhat = U^T*b.
00492      *
00493      * Output: G : GCV function for Tikhonov at alpha
00494      */
00495 {
00496     float t0, G;
00497     float * hold = myvector(m-n);
00498     float * tt= myvector(n);
00499     float * t1= myvector(n);
00500     float * t2= myvector(n);
00501     float * t3= myvector(n);
00502     int i;
00503     
00504     for (i=0; i< m-n; i++)
00505       hold[i] = fabs(bhat[n])*fabs(bhat[n]);
00506     t0 = sum_vector(hold,m-n);
00507 
00508     for (i = 0; i<n; i++){
00509       tt[i] = (float)1.0 /(s[i]*s[i] +alpha*alpha);
00510       t1[i] = alpha*alpha * tt[i];
00511       t2[i] = fabs(bhat[i] * t1[i]) * fabs(bhat[i] * t1[i]);
00512       t3[i] = t1[i] + ((float)1.0 - omega)* (s[i]*s[i]) * tt[i];
00513     }
00514     G = (float) n*(sum_vector(t2,n) + t0) / ((sum_vector(t3,n) + (float)m - (float)n)* (sum_vector(t3,n) + (float)m - (float)n));
00515 
00516 //printf(" svec = %f, to = %f, Num = %f; Dem = %f \n", sum_vector(t2,n), t0, (float) n*(sum_vector(t2,n) + t0), ((sum_vector(t3,n) + (float)m - (float)n)* (sum_vector(t3,n) + (float)m - (float)n)));
00517 
00518     free_vector(hold);
00519     free_vector(tt);
00520     free_vector(t1);
00521     free_vector(t2);
00522     free_vector(t3);
00523     
00524     return G;
00525 }
00526 
00527 /* --------GENERATE_NOISE CODE--------------------------------------------*/
00528 
00529 void generate_noise(float *E, int m, int n)
00530    /*
00531     * Generate random noise.
00532     *
00533     * Input:  
00534     *       m,n : Dimensions of E (needs to be a vector!)
00535     *
00536     * Output: E : Array containing random noise, with norm = 1.
00537     */
00538 {
00539     float **U = matrix(m,12);
00540     float nm, a = -6;
00541     int i, sz = 1, j;
00542     float *o = myvector(m);
00543 
00544     fill_ones(o,m);
00545        
00546     fill_rand(U, m, 12);
00547     row_sum(U, m, 12, E);
00548     
00549     for(j=0; j<m; j++)
00550       E[j] = E[j] + a*o[j];
00551 
00552     //daxpy_(&m, &a, o, &sz, E, &sz);    
00553     //norm = dnrm2_(&m, E, &sz);
00554     nm = norm(E, m);
00555 
00556     for(i = 0; i < m; i++)
00557       E[i] = E[i] / nm;
00558     
00559     free_matrix(U);
00560     free_vector(o);
00561 
00562 }
00563 
00564 
00565 // ======================= imported from matrix.cpp =================
00566 
00567 /* matrix.c contains all the functions that pertain to matrices and vectors
00568  *
00569  * The following are MATRIX functions:
00570  *
00571  *      float **matrix(int m, int n);
00572  *  float **submatrix(float **A, int m, int n, int newm, int newn);
00573  *  void free_matrix(float **A);
00574  *  void fill_rand(float **A, int m, int n);
00575  *  void print_matrix(float **A, int m, int n);
00576  *  void transpose_matrix(float **A, int m, int n);
00577  *  void row_sum(float **A, int m, int n, float *e);
00578  *  void copy_matrix(float **A, float **B, int m, int n);
00579  *
00580  *  The following are VECTOR functions:
00581  *
00582  *  float *myvector(int n);
00583  *  float *subvector(float *x, int n, int newn);
00584  *  void free_vector(float * x);
00585  *  void print_vector(float *x, int n);
00586  *  void fill_ones(float *x, int n);
00587  *  float sum_vector(float *x, int n);
00588  *  void copy_vector(float *x, float *y, int n);
00589  *  void fillveczeros(float *x, int n);
00590  *  float norm(float *x, int n);
00591  *  void alphaxpy(float alpha, float *x, float *y, int n)
00592  *
00593  *  The following are MATRIX-VECTOR functions:
00594  *
00595  *  void matvec_mult(float **A, float *x, int m, int n, float *b);
00596  *  void matvec_multT(float **A, float *x, int m, int n, float *b);
00597  *  float * get_col(float ** A, int m, int k);
00598  *  float * get_row(float ** A, int n, int k);
00599  *  void put_col(float ** A, int m, float * v, int k);
00600  */
00601 
00602 
00603 /* -------------MATRIX FUNCTIONS ------------------------------------*/
00604 
00605 /* Allocates a float matrix A of size m x n with subscript range A[0..m-1][0..n-1] */
00606 float **matrix(int m, int n)
00607 {
00608   int i;
00609   float **A=(float **) calloc(m, sizeof(float*));
00610   if (A == NULL)
00611   {
00612     fprintf(stderr, "Allocation failure in matrix()");
00613     exit (1);
00614   }
00615   A[0]=(float *) calloc(m * n, sizeof(float));
00616   if (A[0] == NULL)
00617   {
00618     fprintf(stderr, "Allocation failure in matrix()");
00619     exit (1);
00620   }
00621   for(i=1; i < m; i++)
00622     A[i]=A[0] + i * n;
00623   return A;
00624 }
00625 
00626 /* Points a submatrix B[0...(newm-1)][0... (newn-1)] to A[0..(newm-1)][0..(newn-1)] */
00627 float **submatrix(float **A, int m, int n, int newm, int newn)
00628 {
00629   int i,j;
00630   float **B = matrix(newm, newn);
00631   if (newm > m || newn > n){
00632     printf("error in submatrix\n");
00633 
00634     return B;
00635   }
00636   for(i=0;i<newm;i++)
00637     for(j=0;j<newn;j++)
00638       B[i][j]=A[i][j];
00639  
00640   return B;
00641 }
00642 
00643 /* Frees a matrix allocated by matrix() */
00644 void free_matrix(float **A)
00645 {
00646   free(A[0]);
00647   free(A);
00648 }
00649 
00650 /* Fills a matrix A[0..m-1][0..n-1] with random numbers. */
00651 void fill_rand(float **A, int m, int n)
00652 {
00653   int i,j;
00654   for(i = 0; i < m; i++) {
00655     for (j = 0; j < n; j++)
00656       A[i][j]=frand();
00657     
00658   }
00659 }
00660 
00661 /* Prints a matrix A[0..m-1][0..n-1] */
00662 void print_matrix(float **A, int m, int n)
00663 {
00664   int i,j;
00665   for (i=0; i< m; i++)
00666   {
00667     for (j=0; j< n; j++)
00668       printf ("%.15f   ", A[i][j]);
00669       putchar ('\n');
00670   }
00671 }
00672 
00673 /* Does a transpose of matrix A and stores it in matrix A*/
00674 void transpose_matrix(float **A, int m, int n)
00675 {
00676   float save;
00677   int i,j;
00678    for (i = 0; i < m; i++){
00679      for (j = i+1; j < n; j++)
00680       {   
00681         save = A[i][j];
00682         A[i][j] = A[j][i];
00683         A[j][i] = save;
00684       }
00685    }
00686 }
00687 
00688 
00689 /* Rearranges matrix A to be in column-major form and stores it in matrix B*/
00690 void column_orient(float **A, int m, int n, float ** B)
00691 {
00692   int i,j;
00693    for (j = 0; j < n; j++){
00694      for (i = 0; i < m; i++)
00695       {   
00696         B[j][i] = A[i][j];
00697       }
00698    }
00699 }
00700 
00701 /* Computes the row sums of matrix a and stores it in vector e*/
00702 void row_sum(float **A, int m, int n, float *e)
00703 {
00704   int i,j;
00705    for (i = 0; i < m; i++){
00706      for (j = 0; j < n; j++)
00707       {   
00708         e[i] = e[i]+ A[i][j];
00709       }
00710    }
00711 }
00712 
00713 /* Copies matrix A into matrix B. */
00714 void copy_matrix(float **A, float **B, int m, int n)
00715 {
00716   int i,j;
00717   for(i = 0; i < m; i++) {
00718     for (j = 0; j < n; j++)
00719       B[i][j]=A[i][j];
00720     
00721   }
00722 }
00723 
00724 
00725 /* -------VECTOR FUNCTIONS ---------------------------------------*/
00726 
00727 /* Allocates a float vector with subscript range x[0..n-1] */
00728 float *myvector(int n)
00729 {
00730   float *x=(float *) calloc(n, sizeof(float));
00731   if (x == NULL)
00732   {
00733     fprintf(stderr, "Allocation failure in matrix()");
00734     exit (1);
00735   }
00736   return x;
00737 }
00738 
00739 /* Points a subvector y[0...(newn-1)] to x[0..(newn-1)] */
00740 float *subvector(float *x, int n, int newn)
00741 {
00742   int i;
00743   float *y = myvector(newn);
00744   if (newn > n){
00745     printf("error in subvector");
00746     return y;
00747   }
00748   for(i=0;i<newn;i++)
00749       y[i]=x[i];
00750   
00751   return y;
00752 }
00753 
00754 /* Frees a vector allocated by vector() */
00755 void free_vector(float *x)
00756 {
00757   free(x);
00758 }
00759 
00760 /* Prints a vector x[0..n-1] */
00761 void print_vector(float *x, int n)
00762 {
00763   int i;
00764   for (i=0; i< n; i++)
00765   {
00766       printf ("%.15f  \n ", x[i]);
00767   }
00768 }
00769 
00770 /* Fills a vector with all ones*/
00771 void fill_ones(float *x, int m)
00772 {
00773   int i;
00774    for (i = 0; i < m; i++){
00775      x[i] = (float) 1.0;
00776    }
00777 }
00778 
00779 /* sums a vector*/
00780 float sum_vector(float *x, int n)
00781 {
00782   int i;
00783   float sum = 0.0;
00784    for (i = 0; i < n; i++){
00785      sum = sum + x[i];
00786    }
00787   return sum;
00788 }
00789 
00790 /* Copies vector x into vector y. */
00791 void copy_vector(float *x, float *y, int n)
00792 {
00793   int i;
00794   for(i = 0; i < n; i++) 
00795       y[i]=x[i];    
00796 }
00797 
00798 /* Fills a vector x with all zeros */
00799 void fillveczeros(float *x, int n)
00800 {
00801  int i;
00802  for (i=0; i<n; i++)
00803    x[i] = 0.0;
00804 }
00805 
00806 /* Computes the Euclidean norm of vector x */
00807 float norm(float *x, int n)
00808 {
00809   float x_norm = 0.0;
00810   int i;
00811   for (i=0; i<n; i++)
00812         x_norm = x_norm + x[i]*x[i];
00813 
00814   x_norm = sqrt(x_norm);
00815 
00816 //x_norm = cblas_dnrm2(n, x, 1);
00817 
00818   return x_norm;
00819 }
00820 
00821 /* Computes a daxpy: y = alpha*x + y  */
00822 void alphaxpy(float alpha, float *x, float *y, int n)
00823 { 
00824  int i;
00825  for (i=0; i<n; i++)
00826    y[i] = alpha*x[i] + y[i];
00827 }
00828 
00829 /* ------MATRIX-VECTOR FUNCTIONS ------------------------------------*/
00830 
00831 /* Does a matrix vector multiply: b = A*x */
00832 void matvec_mult(float **A, float *x, int m, int n, float *b)
00833 {
00834 int i,j;
00835 for(i=0; i< m; i++){
00836         b[i] = 0.0;
00837         for(j=0; j<n; j++){
00838                 b[i] = b[i] + A[i][j]*x[j];
00839 }
00840 }
00841 
00842 //cblas_dgemv(CblasRowMajor, CblasNoTrans, m, n, 1.0, *A, n, x, 1, 0.0, b, 1);
00843 }
00844 
00845 /* Does a matrix-TRANSPOSE vector multiply: b = A'*x */
00846 void matvec_multT(float **A, float *x, int m, int n, float *b)
00847 {
00848         int i,j;
00849 for(i=0; i< n; i++){
00850         b[i] = 0.0;
00851         for(j=0; j<m; j++){
00852                 b[i] = b[i] + A[j][i]*x[j];
00853 }
00854 }
00855 
00856 //cblas_dgemv(CblasRowMajor, CblasTrans, m, n, 1.0, *A, m, x, 1, 0.0, b, 1);
00857 }
00858 
00859 /* The following function gets the kth column of A and 
00860  * copies it into the vector ak, which is returned.
00861  */
00862 void get_col(float ** A, int m, int k, float * ak)
00863 {
00864   int i;
00865   for (i = 0; i < m; i++)
00866     ak[i] = A[i][k];
00867 
00868 }
00869 
00870 /* The following function gets the kth row of A and 
00871  * copies it into the vector ak, which is returned.
00872  */
00873 void get_row(float ** A, int n, int k, float * ak)
00874 {
00875   int i;
00876   for (i = 0; i < n; i++)
00877     ak[i] = A[k][i];
00878 
00879 }
00880 
00881 /* The following function puts vector v into the kth column of A*/
00882 void put_col(float ** A, int m, float * v, int k)
00883 {
00884   int i;
00885   for (i = 0; i < m; i++)
00886     A[i][k] = v[i];
00887 
00888 }
00889