IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 37632


Ignore:
Timestamp:
Nov 19, 2014, 6:46:03 AM (12 years ago)
Author:
eugene
Message:

ill-conditioned matrices in SH fitting...

Location:
branches/eam_branches/ipp-20140904/Ohana/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/gaussj.c

    r36084 r37632  
    11# include <ohana.h>
    2 # define GROWTHTEST 0
     2# define GROWTHTEST 1
    33# define MAX_RANGE 1.0e9
    44
    … …  
    4848    // search for the next pivot
    4949    for (row = 0; row < N; row++) {
    50       if (!finite(A[row][diag])) goto escape;
     50      if (!finite(A[row][diag])) {
     51        fprintf (stderr, "infinity\n");
     52        goto escape;
     53      }
    5154
    5255      // if we have already operated on this row (pivot[row] is true), skip it
    … …  
    6568# if (GROWTHTEST)
    6669    fprintf (stderr, "maxcol: %d\n", maxcol);
    67     fprintf (stderr, "full A matrix:\n");
    68     for (row = 0; row < N; row++) {
    69         for (col = 0; col < N; col++) {
    70             fprintf (stderr, "%10.3e ", A[row][col]);
    71         }
    72         fprintf (stderr, "\n");
    73     }
    74     fprintf (stderr, "\n");
     70    // fprintf (stderr, "full A matrix:\n");
     71    // for (row = 0; row < N; row++) {
     72    //   for (col = 0; col < N; col++) {
     73    //  fprintf (stderr, "%10.3e ", A[row][col]);
     74    //   }
     75    //   fprintf (stderr, "\n");
     76    // }
     77    // fprintf (stderr, "\n");
    7578# endif
    7679
    7780    // if pivot[maxcol] is set, we have already done this row: this implies a singular matrix
    78     if (pivot[maxcol]) goto escape;
     81    if (pivot[maxcol]) {
     82      fprintf (stderr, "singular\n");
     83      goto escape;
     84    }
    7985    pivot[maxcol] = TRUE;
    8086
    … …  
    8692    rowIndex[diag] = maxrow;
    8793    colIndex[diag] = maxcol;
    88     if (A[maxcol][maxcol] == 0.0) goto escape;
     94    if (A[maxcol][maxcol] == 0.0) {
     95      fprintf (stderr, "zero pivot\n");
     96      goto escape;
     97    }
     98
    8999    // XXX Kahan replaces the 0.0 pivot with epsilon*(largest element in column) + underFlow
    90100
    … …  
    102112#   if (GROWTHTEST)
    103113    fprintf (stderr, "column: %d, maxval : %f, growth: %e, epsilon: %e\n", maxcol, tmpval, growth, epsilon);
    104     fprintf (stderr, "A diagonal: ");
    105     for (col = 0; col < N; col++) fprintf (stderr, "%f ", A[col][col]);
    106     fprintf (stderr, "\n");
    107 # endif
    108 
    109     if (fabs(growth) > MAX_RANGE) goto escape;
     114    // fprintf (stderr, "A diagonal: ");
     115    // for (col = 0; col < N; col++) fprintf (stderr, "%f ", A[col][col]);
     116    // fprintf (stderr, "\n");
     117# endif
     118
     119    if (fabs(growth) > MAX_RANGE) {
     120      fprintf (stderr, "max range\n");
     121      goto escape;
     122    }
    110123
    111124    /* adjust the elements above the pivot */
    … …  
    120133
    121134# if (GROWTHTEST)
    122     fprintf (stderr, "final A matrix:\n");
    123     for (row = 0; row < N; row++) {
    124         for (col = 0; col < N; col++) {
    125             fprintf (stderr, "%10.3e ", A[row][col]);
    126         }
    127         fprintf (stderr, "\n");
    128     }
    129     fprintf (stderr, "\n");
     135  fprintf (stderr, "final A matrix:\n");
     136  // for (row = 0; row < N; row++) {
     137  //   for (col = 0; col < N; col++) {
     138  //     fprintf (stderr, "%10.3e ", A[row][col]);
     139  //   }
     140  //   fprintf (stderr, "\n");
     141  // }
     142  // fprintf (stderr, "\n");
    130143# endif
    131144

  • branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/ohana_allocate.c

    r37547 r37632  
    181181
    182182  ref = (Memblock *) in - 1;
     183
     184  // fprintf (stderr, "  file: %s, line: %d, func: %s, size: %zd, addr: %zx\n",
     185  // ref->file, ref->line, ref->func, ref->size, (size_t) ref);
    183186
    184187  pthread_mutex_lock(&memBlockListMutex);

  • branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/shfit.c

    r37630 r37632  
    3434  ResetVector (Mvec,  OPIHI_FLT, terms->Nterms);
    3535
    36   double **Ar, **br, **Ai, **bi;
    37   ALLOCATE (Ar, double *, terms->Nterms);
    38   ALLOCATE (Ai, double *, terms->Nterms);
    39   ALLOCATE (br, double *, terms->Nterms);
    40   ALLOCATE (bi, double *, terms->Nterms);
     36  // we only fit the linearlly independent terms: Re(m >= 0), Im(m > 0)
     37  int Nre = 0;
     38  int Nim = 0;
    4139  for (i = 0; i < terms->Nterms; i++) {
    42     ALLOCATE_ZERO (Ar[i], double, terms->Nterms);
    43     ALLOCATE_ZERO (Ai[i], double, terms->Nterms);
    44     ALLOCATE_ZERO (br[i], double, 1);
    45     ALLOCATE_ZERO (bi[i], double, 1);
     40    if (terms->m[i] >= 0) Nre ++;
     41    if (terms->m[i] >  0) Nim ++;
     42  }
     43  int *Jre = NULL;
     44  int *Jim = NULL;
     45  ALLOCATE (Jre, int, Nre);
     46  ALLOCATE (Jim, int, Nim);
     47
     48  Nre = Nim = 0;
     49  for (i = 0; i < terms->Nterms; i++) {
     50    if (terms->m[i] >= 0) {
     51      Jre[Nre] = i;
     52      Nre ++;
     53    }
     54    if (terms->m[i] >  0) {
     55      Jim[Nim] = i;
     56      Nim ++;
     57    }
     58  }
     59
     60  double **Are, **bre, **Aim, **bim;
     61  ALLOCATE (Are, double *, Nre);
     62  ALLOCATE (bre, double *, Nre);
     63  ALLOCATE (Aim, double *, Nim);
     64  ALLOCATE (bim, double *, Nim);
     65  for (i = 0; i < Nre; i++) {
     66    ALLOCATE_ZERO (Are[i], double, Nre);
     67    ALLOCATE_ZERO (bre[i], double, 1);
     68  }
     69  for (i = 0; i < Nim; i++) {
     70    ALLOCATE_ZERO (Aim[i], double, Nim);
     71    ALLOCATE_ZERO (bim[i], double, 1);
    4672  }
    4773
    … …  
    5379    double Fv = Fvec->elements.Flt[i];
    5480
    55     for (j = 0; j < terms->Nterms; j++) {
    56       br[j][0] += Fv * terms->Fr[j];
    57       bi[j][0] += Fv * terms->Fi[j];
     81    for (j = 0; j < Nre; j++) {
     82      int jre = Jre[j];
     83      bre[j][0] += Fv * terms->Fr[jre];
     84    }
     85    for (j = 0; j < Nim; j++) {
     86      int jim = Jim[j];
     87      bim[j][0] += Fv * terms->Fi[jim];
    5888    }
    5989
    60     for (j = 0; j < terms->Nterms; j++) {
    61       for (k = j; k < terms->Nterms; k++) {
    62         Ar[j][k] += terms->Fr[j] * terms->Fr[k];
    63         Ai[j][k] += terms->Fi[j] * terms->Fi[k];
    64       }
     90    for (j = 0; j < Nre; j++) {
     91      int jre = Jre[j];
     92      for (k = j; k < Nre; k++) {
     93        int kre = Jre[k];
     94        Are[j][k] += terms->Fr[jre] * terms->Fr[kre];
     95      }
    6596    }
    66   }
    6797
    68   for (j = 1; j < terms->Nterms; j++) {
    69     for (k = 0; k < j; k++) {
    70       Ar[j][k] = Ar[k][j];
    71       Ai[j][k] = Ai[k][j];
    72     }   
    73   }
    74 
    75   for (j = 0; j < terms->Nterms; j++) {
    76     if (terms->m[j] == 0) {
    77       Ai[j][j] = 1.0; // Im(Y_l0) == 0.0
    78       bi[j][0] = 0.0;
    79     } else {
    80       if (terms->m[j] < 0) {
    81         for (k = 0; k < terms->Nterms; k++) {
    82           Ar[j][k] = 0.0;
    83           Ar[k][j] = 0.0;
    84           Ai[j][k] = 0.0;
    85           Ai[k][j] = 0.0;
    86         }
    87         Ar[j][j] = 1.0;
    88         Ai[j][j] = 1.0;
    89         br[j][0] = 0.0;
    90         bi[j][0] = 0.0;
    91         // Re(Y_lm) = Re(Y_l(-m)) [m even], Re(Y_lm) = Re(-Y_l(-m)) [m odd]
     98    for (j = 0; j < Nim; j++) {
     99      int jim = Jim[j];
     100      for (k = j; k < Nim; k++) {
     101        int kim = Jim[k];
     102        Aim[j][k] += terms->Fi[jim] * terms->Fi[kim];
    92103      }
    93104    }
    94105  }
    95106
    96   if (0) {
    97     fprintf (stderr, "--- Ar --- : br \n");
    98     for (i = 0; i < terms->Nterms; i++) {
    99       for (j = 0; j < terms->Nterms; j++) {
    100         fprintf (stderr, "%10.6f ", Ar[i][j]);
    101       }
    102       fprintf (stderr, " : %10.6f\n", br[i][0]);
    103     }
    104 
    105     fprintf (stderr, "--- Ai --- : bi \n");
    106     for (i = 0; i < terms->Nterms; i++) {
    107       for (j = 0; j < terms->Nterms; j++) {
    108         fprintf (stderr, "%10.6f ", Ai[i][j]);
    109       }
    110       fprintf (stderr, " : %10.6f\n", bi[i][0]);
    111     }
     107  for (j = 1; j < Nre; j++) {
     108    for (k = 0; k < j; k++) {
     109      Are[j][k] = Are[k][j];
     110    }   
     111  }
     112  for (j = 1; j < Nim; j++) {
     113    for (k = 0; k < j; k++) {
     114      Aim[j][k] = Aim[k][j];
     115    }   
    112116  }
    113117
    114   fprintf (stderr, "--- Ar --- : br \n");
    115   for (j = 0; j < terms->Nterms; j++) {
    116     fprintf (stderr, "%10.6f ", Ar[j][j]);
    117     fprintf (stderr, " : %10.6f\n", br[j][0]);
     118  fprintf (stderr, "--- Are --- : bre \n");
     119  for (j = 0; j < Nre; j++) {
     120    fprintf (stderr, "%10.6f ", Are[j][j]);
     121    fprintf (stderr, " : %10.6f\n", bre[j][0]);
    118122  }
    119123
    120   fprintf (stderr, "--- Ai --- : bi \n");
    121   for (i = 0; i < terms->Nterms; i++) {
    122     fprintf (stderr, "%10.6f : ", Ai[i][i]);
    123     fprintf (stderr, " : %10.6f\n", bi[i][0]);
     124  fprintf (stderr, "--- Aim --- : bim \n");
     125  for (i = 0; i < Nim; i++) {
     126    fprintf (stderr, "%10.6f : ", Aim[i][i]);
     127    fprintf (stderr, " : %10.6f\n", bim[i][0]);
    124128  }
    125129
    126   if (!dgaussjordan (Ar, br, terms->Nterms, 1)) {
     130  if (!dgaussjordan (Are, bre, Nre, 1)) {
    127131    gprint (GP_ERR, "failed to fit data : ill-conditioned matrix\n");
    128132    return FALSE;
    129133  }
    130   if (!dgaussjordan (Ai, bi, terms->Nterms, 1)) {
     134  if (!dgaussjordan (Aim, bim, Nim, 1)) {
    131135    gprint (GP_ERR, "failed to fit data : ill-conditioned matrix\n");
    132136    return FALSE;
    133137  }
    134138
    135   for (j = 0; j < terms->Nterms; j++) {
    136     Fr[j] = br[j][0] / (float) Frvec[0].Nelements;
    137     Fi[j] = bi[j][0] / (float) Fivec[0].Nelements;
     139  // XXX EAM : why is this factor needed?
     140  for (j = 0; j < terms->Nterms; j++) {
     141    Fr[j] = 0.0;
     142    Fi[i] = 0.0;
     143  }
     144  for (j = 0; j < Nre; j++) {
     145    int jre = Jre[j];
     146    Fr[jre] = bre[j][0] * terms->Nterms / (float) Frvec[0].Nelements;
     147  }
     148  for (j = 0; j < Nim; j++) {
     149    int jim = Jim[j];
     150    Fi[jim] = bim[j][0] * terms->Nterms / (float) Frvec[0].Nelements;
    138151  }
    139152
    … …  
    142155    Lvec[0].elements.Flt[j]  = terms->l[j];
    143156    Mvec[0].elements.Flt[j]  = terms->m[j];
    144     Frvec[0].elements.Flt[j] = 4*M_PI*Fr[j];
    145     Fivec[0].elements.Flt[j] = 4*M_PI*Fi[j];
     157    Frvec[0].elements.Flt[j] = Fr[j];
     158    Fivec[0].elements.Flt[j] = Fi[j];
    146159  }
    147160
Note: See TracChangeset for help on using the changeset viewer.