cospmls.cpp

普林斯顿开发的快速球面调和变换算法

  /* note that the number of non-zero entries is the same
     as in the non-transposed case */

  trans_tableptr = result;
  
  /* now traverse the cos_pml_table , loading appropriate values
     into the rows of transposed array */

  for (row = 0; row < bw; row++)
    {
      /* if m odd, no need to do last row - all zeroes */
      if (row == (bw-1))
	  {
	    if ((m % 2) == 1)
	      return;
	  }

      /* get the rowsize for the transposed array */
      rowsize = Transpose_RowSize(row, m, bw);

      /* compute the starting point for values in cos_pml_table */
      if (row <= m)
	{
	  if ((row % 2) == 0)
	    tableptr = cos_pml_table + (row/2);
	  else
	    tableptr = cos_pml_table + (m/2) + 1 + (row/2);
	}
      else
	{
	  /* if row > m, then the highest degree coefficient
	     of P(m,row) should be the first coefficient loaded
	     into the transposed array, so figure out where
	     this point is.
	  */
	  offset = 0;
	  if ((m%2) == 0)
	    {
	      for (i=m; i<=row; i++)
		offset += RowSize(m, i);
	    }
	  else
	    {
	      for (i=m;i<=row+1;i++)
		offset += RowSize(m, i);
	    }
	  /* now we are pointing one element too far, so decrement */
	  offset--;

	  tableptr = cos_pml_table + offset;
	}

      /* stride is how far we need to jump between
	 values in cos_pml_table, i.e., to traverse the columns of the
	 cos_pml_table.  Need to set initial value.  Stride always
	 increases by 2 after that 
      */
      if (row <= m)
	stride = m + 2 - (m % 2) + (row % 2);
      else
	stride = row + 2;

      /* now load up this row of the transposed table */
      costable_offset = 0;
      for (i=0; i < rowsize; i++)
	{
	  trans_tableptr[i] = tableptr[costable_offset];
	  costable_offset += stride;
	  stride += 2;

	} /* closes i loop */

      trans_tableptr += rowsize;

    } /* closes row loop */


}

double** cospmls::Spharmonic_Pml_Table(int bw, double *resultspace, double *workspace)
{
  int i;
  double **spharmonic_pml_table;

  /* allocate an array of double pointers */
  spharmonic_pml_table = (double **) malloc(sizeof(double *) * bw);

  /* traverse the array, assigning a location in the resultspace
     to each pointer */

  spharmonic_pml_table[0] = resultspace;

  for (i=1; i<bw; i++)
    {
      spharmonic_pml_table[i] = spharmonic_pml_table[i-1] +
	                        this->TableSize(i-1,bw);
    }
  
  /* now load up the array with CosPml and CosGml values */
  for (i=0; i<bw; i++)
    {
      this->CosPmlTableGen(bw, i, spharmonic_pml_table[i], workspace);
    }

  /* that's it */

  return spharmonic_pml_table;
}

double** cospmls::Reduced_Spharmonic_Pml_Table(int bw, int m, double *resultspace, double *workspace)
{
	 int i;
  double **spharmonic_pml_table;

  /* allocate an array of double pointers */
  spharmonic_pml_table = (double **) malloc(sizeof(double *) * (m + 1));

  /* traverse the array, assigning a location in the resultspace
     to each pointer */

  spharmonic_pml_table[0] = resultspace;

  for (i=1; i<m; i++)
    {
      spharmonic_pml_table[i] = spharmonic_pml_table[i-1] +
	                        TableSize(i-1,bw);
    }
  
  /* now load up the array with CosPml and CosGml values */
  for (i=0; i<m; i++)
    {
      CosPmlTableGen(bw, i, spharmonic_pml_table[i], workspace);
    }

  /* that's it */

  return spharmonic_pml_table;



}

double** cospmls::Transpose_Spharmonic_Pml_Table(double **spharmonic_pml_table, int bw, double *resultspace, double *workspace)
{
  int i;
  double **transpose_spharmonic_pml_table;

  /* allocate an array of double pointers */
  transpose_spharmonic_pml_table = (double **) malloc(sizeof(double *) * bw);

  /* now need to load up the transpose_spharmonic_pml_table by transposing
     the tables in the spharmonic_pml_table */

  transpose_spharmonic_pml_table[0] = resultspace;

  for (i=0; i<bw; i++)
    {
      Transpose_CosPmlTableGen(bw, 
			       i, 
			       spharmonic_pml_table[i],
			       transpose_spharmonic_pml_table[i]);

      if (i != (bw-1))
	  {
	    transpose_spharmonic_pml_table[i+1] =
	      transpose_spharmonic_pml_table[i] + TableSize(i, bw);
	  }
    }

  return transpose_spharmonic_pml_table;
	

}

double** cospmls::Reduced_Transpose_Spharmonic_Pml_Table(double **spharmonic_pml_table, int bw, int m, double *resultspace, double *workspace)
{
	 int i;
  double **transpose_spharmonic_pml_table;

  /* allocate an array of double pointers */
  transpose_spharmonic_pml_table = (double **) malloc(sizeof(double *) * m);

  /* now need to load up the transpose_spharmonic_pml_table by transposing
     the tables in the spharmonic_pml_table */

  transpose_spharmonic_pml_table[0] = resultspace;

  for (i=0; i<m; i++)
    {
      Transpose_CosPmlTableGen(bw, 
			       i, 
			       spharmonic_pml_table[i],
			       transpose_spharmonic_pml_table[i]);

      if (i != (bw-1))
	  {
	    transpose_spharmonic_pml_table[i+1] =
	      transpose_spharmonic_pml_table[i] + TableSize(i, bw);
	  }
    }

  return transpose_spharmonic_pml_table;


}
void cospmls::PmlTableGen(int bw, int m, double *storeplm, double *workspace)
{
  double *prev, *prevprev, *temp1, *temp2, *temp3, *temp4, *x_i, *eval_args;
  int i;
  /* double *storeplm_ptr; */
  
  prevprev = workspace;
  prev = prevprev + (2*bw);
  temp1 = prev + (2*bw);
  temp2 = temp1 + (2*bw);
  temp3 = temp2 + (2*bw);
  temp4 = temp3 + (2*bw);
  x_i = temp4 + (2*bw);
  eval_args = x_i + (2*bw);
  
  /* storeplm_ptr = storeplm; */

  /* get the evaluation nodes */
  prim.EvalPts(2*bw,x_i);
  prim.ArcCosEvalPts(2*bw,eval_args);
  
  /* set initial values of first two Pmls */
  for (i=0; i<2*bw; i++) 
    prevprev[i] = 0.0;
  if (m == 0)
    for (i=0; i<2*bw; i++)
      prev[i] = 0.5;
  else 
    prim.Pmm_L2(m, eval_args, 2*bw, prev);

  memcpy(storeplm, prev, (size_t) sizeof(double) * 2 * bw);

  for(i = 0; i < bw - m - 1; i++)
    {
      prim.vec_mul(prim.L2_cn(m,m+i),prevprev,temp1,2*bw);
      prim.vec_pt_mul(prev, x_i, temp2, 2*bw);
      prim.vec_mul(prim.L2_an(m,m+i), temp2, temp3, 2*bw);
      prim.vec_add(temp3, temp1, temp4, 2*bw); /* temp4 now contains P(m,m+i+1) */
      
      storeplm += (2 * bw);
      memcpy(storeplm, temp4, (size_t) sizeof(double) * 2 * bw);
      memcpy(prevprev, prev, (size_t) sizeof(double) * 2 * bw);
      memcpy(prev, temp4, (size_t) sizeof(double) * 2 * bw);
    }

  /*  storeplm = storeplm_ptr; */



}

int cospmls::Reduced_Naive_TableSize(int bw, int m)
{

	 
  int i, sum;

  sum = 0;
  
  for (i=m; i<bw; i++)
      sum += ( 2 * bw * (bw - i));

  return sum;


}

double** cospmls::SemiNaive_Naive_Pml_Table(int bw, int m, double *resultspace, double *workspace)
{
	 int i;
  double **seminaive_naive_table;
  int  lastspace;

  seminaive_naive_table = (double **) malloc(sizeof(double) * (bw+1));

  seminaive_naive_table[0] = resultspace;

  
  for (i=1; i<m; i++)
    {
      seminaive_naive_table[i] = seminaive_naive_table[i - 1] +
	                        TableSize(i-1,bw);
    }

  if( m == 0)
    {
      lastspace = 0;
      for (i=m+1; i<bw; i++)
	{ 
	  seminaive_naive_table[i] = seminaive_naive_table[i - 1] +
	    (2 * bw * (bw - (i - 1)));
	}
    }
  else
    {
      lastspace = TableSize(m-1,bw);
      seminaive_naive_table[m] = seminaive_naive_table[m-1] +
	lastspace;
      for (i=m+1; i<bw; i++)
	{ 
	  seminaive_naive_table[i] = seminaive_naive_table[i - 1] +
	    (2 * bw * (bw - (i - 1)));
	}
    }

  /* now load up the array with CosPml and CosGml values */
  for (i=0; i<m; i++)
    {
      CosPmlTableGen(bw, i, seminaive_naive_table[i], workspace);
    }

  /* now load up pml values */
  for(i=m; i<bw; i++)
    {
      PmlTableGen(bw, i, seminaive_naive_table[i], workspace);
    }

  /* that's it */

  return seminaive_naive_table;


}

double** cospmls::Transpose_SemiNaive_Naive_Pml_Table(double **seminaive_naive_pml_table, int bw, int m, double *resultspace, double *workspace)
{
	 int i, lastspace;
  double **trans_seminaive_naive_pml_table;

  /* allocate an array of double pointers */
  trans_seminaive_naive_pml_table = (double **) malloc(sizeof(double *) * (bw+1));

  /* now need to load up the transpose_seminaive_naive_pml_table by transposing
     the tables in the seminiave portion of seminaive_naive_pml_table */

  trans_seminaive_naive_pml_table[0] = resultspace;

  
  for (i=1; i<m; i++)
    {
      trans_seminaive_naive_pml_table[i] =
	trans_seminaive_naive_pml_table[i - 1] +
	TableSize(i-1,bw);
    }

  if( m == 0)
    {
      lastspace = 0;
      for (i=m+1; i<bw; i++)
	{ 
	  trans_seminaive_naive_pml_table[i] =
	    trans_seminaive_naive_pml_table[i - 1] +
	    (2 * bw * (bw - (i - 1)));
	}
    }
  else
    {
      lastspace = TableSize(m-1,bw);
      trans_seminaive_naive_pml_table[m] =
	trans_seminaive_naive_pml_table[m-1] +
	lastspace;

      for (i=m+1; i<bw; i++)
	{ 
	  trans_seminaive_naive_pml_table[i] =
	    trans_seminaive_naive_pml_table[i - 1] +
	    (2 * bw * (bw - (i - 1)));
	}
    }

  for (i=0; i<m; i++)
    {
      Transpose_CosPmlTableGen(bw, 
			       i, 
			       seminaive_naive_pml_table[i],
			       trans_seminaive_naive_pml_table[i]);

      if (i != (bw-1))
	  {
	    trans_seminaive_naive_pml_table[i+1] =
	      trans_seminaive_naive_pml_table[i] + TableSize(i, bw);
	  }
    }

  /* now load up pml values */
  for(i=m; i<bw; i++)
    {
      PmlTableGen(bw, i, trans_seminaive_naive_pml_table[i], workspace);
    }

  return trans_seminaive_naive_pml_table;

}