cospmls.c

基于java的3d开发库。对坐java3d的朋友有很大的帮助。

  */
  
  double *trans_tableptr, *tableptr;
  int i, row, rowsize, stride, offset, costable_offset;

  /* 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 */

}
/************************************************************************/
/* This is a function that returns all of the (cosine transforms of)
   Pmls and Gmls necessary
   to do a full spherical harmonic transform, i.e., it calls
   CosPmlTableGen for each value of m less than bw, returning a
   table of tables ( a pointer of type (double **), which points
   to an array of size m, each containing a (double *) pointer
   to a set of cospml or cosgml values, which are the (decimated)
   cosine series representations of Pml (even m) or Gml (odd m)
   functions.  See CosPmlTableGen for further clarification.

   Inputs - the bandwidth bw of the problem
   resultspace - need to allocate Spharmonic_TableSize(bw) for storing results
   workspace - needs to be (16 * bw)

   Note that resultspace is necessary and contains the results/values
   so one should be careful about when it is OK to re-use this space.
   workspace, though, does not have any meaning after this function is
   finished executing

*/

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

  /* that's it */

  return spharmonic_pml_table;
}


/************************************************************************/
/*  This is a reduced version of Spharmonic_Pml_Table(), reduced in
    the sense that I'm interested in generating tables only up to
    a given order m. That's because this table will be used by
    the hybrid algorithm's call of the semi-naive algorithm, which
    doesn't need the table for all orders. Hopefully, this'll cut
    down on memory usage.

    See documentation of Spharmonic_Pml_Table() for details on how
    what's going on in this function.

   Inputs - the bandwidth bw of the problem, m the maximum order
            one is interested in
   resultspace - need to allocate Reduced_Spharmonic_TableSize(bw,m) for
            storing results
   workspace - needs to be (16 * bw)

   Note that resultspace is necessary and contains the results/values
   so one should be careful about when it is OK to re-use this space.
   workspace, though, does not have any meaning after this function is
   finished executing

*/

double **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;

}

/************************************************************************/
/* For the inverse semi-naive spharmonic transform, need the "transpose"
   of the spharmonic_pml_table.  Need to be careful because the
   entries in the spharmonic_pml_table have been decimated, i.e.,
   the zeroes have been stripped out.

   Inputs are a spharmonic_pml_table generated by Spharmonic_Pml_Table
   and the bandwidth bw

   Allocates memory for the (double **) result
   also allocates memory

   resultspace - need to allocate Spharmonic_TableSize(bw) for storing results
   workspace - not needed, but argument added to avoid
               confusion wth Spharmonic_Pml_Table

*/

double **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;
}


/************************************************************************/
/* For the inverse semi-naive spharmonic transform, need the "transpose"
   of the spharmonic_pml_table.  Need to be careful because the
   entries in the spharmonic_pml_table have been decimated, i.e.,
   the zeroes have been stripped out.

   Inputs are a spharmonic_pml_table generated by Spharmonic_Pml_Table
   and the bandwidth bw, and the order m which is the last order that
   semi-naive is used when computing a full spherical transform

   Allocates memory for the (double **) result
   also allocates memory

   resultspace - need to allocate Reduced_Spharmonic_TableSize(bw,m)
                 for storing results
   workspace - not needed, but argument added to avoid
               confusion wth Reduced_Spharmonic_Pml_Table

   This is a reduced version of Transpose_Spharmonic_Pml_Table.

*/

double **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;
}

/************************************************************************/
/* generate all of the Pmls for a specified value of m.  

   storeplm points to a double array of size 2 * bw * (bw - m);

   Workspace needs to be
   16 * bw 

   P(m,l,j) respresents the associated Legendre function P_l^m
   evaluated at the j-th Chebyshev point (for the bandwidth bw)
   Cos((2 * j + 1) * PI / (2 * bw)).

   The array is placed in storeplm as follows:

   P(m,m,0)    P(m,m,1)  ... P(m,m,2*bw)
   P(m,m+1,0)  P(m,m+1,1)... P(m,m+1,2*bw)
   P(m,m+2,0)  P(m,m+2,1)... P(m,m+2,2*bw)
   ...
   P(m,bw,0)   P(m,bw,1) ... P(m,bw,2*bw)

   This array will eventually be used by the naive transform algorithm.
   This function will precompute the arrays necessary for the algorithm.
*/

static void 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 */
  EvalPts(2*bw,x_i);
  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 
    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++)
    {
      vec_mul(L2_cn(m,m+i),prevprev,temp1,2*bw);
      vec_pt_mul(prev, x_i, temp2, 2*bw);
      vec_mul(L2_an(m,m+i), temp2, temp3, 2*bw);
      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; */
}



/************************************************************************/
/* Reduced_Naive_TableSize(bw,m) returns an integer value for
   the amount of space necessary to fill out a reduced naive table
   of pmls if interested in using it only for orders m through bw - 1.

*/

int 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;

}

/*************************************************************

  just like Spharmonic_Pml_Table(), except generates a
  table for use with the semi-naive and naive algorithms.

  m is the cutoff order, where to switch from semi-naive to
  naive algorithms

  bw = bandwidth of problem
  m  = where to switch algorithms (order where naive is FIRST done)
  resultspace = where to store results, must be of
                size Reduced_Naive_TableSize(bw, m) +
		     Reduced_SpharmonicTableSize(bw, m);

***********************************************************/

double **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;

}


/************************************************************************/
/* For the inverse seminaive_naive transform, need the "transpose"
   of the seminaive_naive_pml_table.  Need to be careful because the
   entries in the seminaive portion have been decimated, i.e.,
   the zeroes have been stripped out.

   Inputs are a seminaive_naive_pml_table generated by SemiNaive_Naive_Pml_Table
   and the bandwidth bw and cutoff order m

   Allocates memory for the (double **) result
   also allocates memory

   resultspace - need to allocate Reduced_Naive_TableSize(bw, m) +
		     Reduced_SpharmonicTableSize(bw, m) for storing results
   workspace - size 16 * bw 

*/

double **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;
}