00001 #ifndef F77_MATRIX_H
00002 #define F77_MATRIX_H
00003 
00004 
00005 #include <assert.h>
00006 #include <sys/types.h>
00007 
00008 
00009 
00010 
00011 
00012 
00013 
00014 
00015 
00016 
00017 
00018 
00019 
00020 
00021 
00022 
00023 
00024 
00025 
00026 
00027 
00028 
00029 
00030 
00031 
00032 template <class T>
00033 class FMATRIX {
00034   public:
00035    FMATRIX(size_t dim1, size_t dim2);
00036    FMATRIX(T* cpparr, size_t dim1, size_t dim2);
00037    operator T*();
00038    T& operator()(size_t index1, size_t index2);
00039    ~FMATRIX();
00040   public:
00041    const size_t ndim;  
00042    size_t dim[7];      
00043    T*  cpprep;         
00044    T*  f77rep;         
00045 };
00046 
00047 template <class T>
00048 FMATRIX<T>::FMATRIX(size_t dim1, size_t dim2)
00049    : cpprep(NULL),f77rep(new T[dim1*dim2]),ndim(2)
00050 {
00051    dim[0]=dim1;
00052    dim[1]=dim2;
00053    dim[2]=0;
00054    dim[3]=0;
00055    dim[4]=0;
00056    dim[5]=0;
00057    dim[6]=0;
00058 }
00059 
00060 template <class T>
00061 FMATRIX<T>::FMATRIX(T* cpparr, size_t dim1, size_t dim2)
00062    : cpprep(cpparr),f77rep(new T[dim1*dim2]),ndim(2)
00063 {
00064    dim[0]=dim1;
00065    dim[1]=dim2;
00066    dim[2]=0;
00067    dim[3]=0;
00068    dim[4]=0;
00069    dim[5]=0;
00070    dim[6]=0;
00071 
00072    
00073    size_t index_cpp=0;
00074    size_t index_f77;
00075    for(size_t i=0;i<dim[0];i++) {
00076       for(size_t j=0;j<dim[1];j++) {
00077          index_f77 = j*dim[0] + i;
00078          f77rep[index_f77] = cpprep[index_cpp++];
00079       }
00080    }
00081 }
00082 
00083 template <class T>
00084 FMATRIX<T>::operator T*()
00085 {
00086    
00087    return f77rep;
00088 }
00089 
00090 template <class T>
00091 T&  FMATRIX<T>::operator()(size_t index1, size_t index2)
00092 {
00093    assert(ndim==2);  
00094 
00095    
00096    size_t index_f77 = index2*dim[0] + index1;
00097 
00098    
00099    return *(f77rep+index_f77);
00100 }
00101 
00102 template <class T>
00103 FMATRIX<T>::~FMATRIX()
00104 {
00105    if(cpprep) {
00106       assert(ndim==2);  
00107 
00108       
00109       size_t index_cpp;
00110       size_t index_f77=0;
00111       for(size_t j=0;j<dim[1];j++) {
00112          for(size_t i=0;i<dim[0];i++) {
00113             index_cpp = i*dim[1] + j;
00114             cpprep[index_cpp] = f77rep[index_f77++];
00115          }
00116       }
00117    }
00118 
00119    
00120    delete[] f77rep;
00121 }
00122 
00123 
00124 #endif