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