transpose.h

The Spectral Toolkit is a C++ spectral transform library written by Rodney James and Chuck Panaccion

//
// spectral toolkit 
// copyright (c) 2005 university corporation for atmospheric research
// licensed under the gnu general public license
//

#ifndef __transpose__
#define __transpose__

#include "spectral.h"

namespace spectral
{

  /// Matrix transpose function template.  Performs blocked transpose of 
  /// submatrix ranging from rows p1 to p1+m1, multiplying output matrix 
  /// by constant c.  
  /// \param a input matrix
  /// \param b output matrix
  /// \param m1 length of row segment 
  /// \param p1 offset of row segment
  /// \param n2 column length
  /// \param c real normalization constant
  
  template <typename T,typename U> 
  void transpose(T **a,T **b,int m1,int p1,int n2,U c)
  {
    int I1,I2;
    int i1,i2;
    const int block=32;
    for(I1=p1;I1<p1+m1;I1+=block)
      for(I2=0;I2<n2;I2+=block)
	for(i1=I1;i1<min(I1+block,p1+m1);i1++)
	  for(i2=I2;i2<min(I2+block,n2);i2++)
	    b[i2][i1]=c*a[i1][i2];
  }
	
  /// Matrix transpose function template.  Performs blocked transpose of 
  /// submatrix ranging from rows p1 to p1+m1, multiplying output matrix 
  /// by constant c.  
  /// \param a input matrix
  /// \param b output matrix
  /// \param m1 length of row segment 
  /// \param p1 offset of row segment
  /// \param n2 column length
  
  template <typename T> 
  void transpose(T **a,T **b,int m1,int p1,int n2)
  {
    int I1,I2;
    int i1,i2;
    const int block=32;
    for(I1=p1;I1<p1+m1;I1+=block)
      for(I2=0;I2<n2;I2+=block)
	for(i1=I1;i1<min(I1+block,p1+m1);i1++)
	  for(i2=I2;i2<min(I2+block,n2);i2++)
	    b[i2][i1]=a[i1][i2];
  }
	

  /// Insitu matrix transpose function template.  Performs insitu transpose of 
  /// m by n matrix, returning n by m matrix.
  /// \param a matrix reference (m by n on input, n by m on output)
  /// \param m row length
  /// \param n column length
  
  template <typename T>
  void transpose(T** &a,int m,int n)
  {
    int i,j;
    T x;
    T *A=a[0];
    int N=m*n-1;
    for(i=1;i<N-1;i++)
      {
	j=n*i%N;
	while(j<i)
	  j=n*j%N;
	if(j!=i) 
	  {
	    x=A[i];
	    A[i]=A[j];
	    A[j]=x;
	  }
      }
    delete [] a;
    a=new T*[n];
    a[0]=A;
    for(int i=1;i<n;i++)
      a[i]=a[i-1]+m;
  }

}

#endif

// Local Variables:
// mode:C++
// End: