home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / listings / v_09_01 / 9n01128a < prev    next >
Text File  |  1990-03-26  |  5KB  |  189 lines

  1. #include <stdlib.h>
  2. #include "cmatrix.hpp"
  3.  
  4. error(char* msg,int index,int s)
  5.    {
  6.    printf("%s%d %d\n",msg,index,s);
  7.    exit(1);
  8.    }
  9.  
  10. void Cvector::check( int index)
  11.    {int loc;
  12.    loc=index-(*this).base;
  13.    if (loc<0)
  14.       error(" Cvector error: index-base<0",index,(*this).base);
  15.    if (loc>= (*this).size)
  16.       error(" Cvector error: index too large",index,(*this).size);
  17.    }
  18.  
  19. void Cmatrix::check( int i,int j)
  20.    {int loc;
  21.    loc=i-(*this).base;
  22.    if (loc<0)
  23.       error(" Cmatrix error: row index-base<0",i,(*this).base);
  24.    if (loc>= (*this).rowkt)
  25.       error(" Cmatrix error: row index too large",i,(*this).rowkt);
  26.  
  27.    loc=j-(*this).base;
  28.    if (loc<0)
  29.       error(" Cmatrix error: column index-base<0",j,(*this).base);
  30.    if (loc>= (*this).colkt)
  31.       error(" Cmatrix error: column index too large",j,(*this).colkt);
  32.  
  33.    }
  34.  
  35.  
  36. void Cmatrix::init(int row, int col, int b)
  37.    {
  38.    if(row<1||col<1)error(" matrix needs at least one row/col ",row,col);
  39.    rowkt=row;colkt=col;base=b;
  40.    m= new complex *[rowkt];
  41.    for(int j=0;j<rowkt;j++)m[j]= new complex[colkt];
  42.    complex zero;complex one(1.,0.);
  43.    //initialize to identity matrix. 
  44.    for(int i=0;i<rowkt;i++){for(j=0;j<colkt;j++)m[i][j]=zero;
  45.       m[i][i]= one;// omit this line for init. to zero matrix
  46.       }
  47.    }
  48.  
  49. Cmatrix::Cmatrix(Cmatrix& a)
  50.    { init(a.rowkt,a.colkt,a.base);
  51.     for(int i=0;i<a.rowkt;i++)for(int j=0;j<a.colkt;j++)m[i][j]=a.m[i][j];
  52.    }
  53.  
  54. Cmatrix::~Cmatrix()
  55.    {for(int i=0;i<rowkt;i++)delete m[i];
  56.    delete m;
  57.    }
  58.  
  59.  
  60. complex& Cvector::element(int i)
  61.    {
  62.    check(i);return head[i-base];
  63.    }
  64.  
  65. void Cvector::setelement( int i, complex& value)
  66.    {
  67.    check(i);
  68.    head[i]=value;
  69.    }
  70. complex& Cmatrix::element(int i,int j)
  71.    {
  72.    check(i,j);return m[i-base][j-base];
  73.    }
  74.  
  75. void Cmatrix::setelement( int i, int j, complex value)
  76.    {
  77.    check(i,j);m[i][j]=value;
  78.    }
  79.  
  80.  
  81. void Cmatrix::operator=(Cmatrix& a)
  82.    {
  83.    if(this==&a)return;
  84.    for(int i=0;i<rowkt;i++)delete m[i];
  85.    delete m;
  86.    init(a.rowkt,a.colkt,a.base);
  87.    for(i=0;i<a.rowkt;i++)for(int j=0;j<a.colkt;j++)m[i][j]=a.m[i][j];
  88.    }
  89.  
  90.  Cvector::Cvector( Cvector& orig) //copy
  91.    {
  92.    size=orig.size;base=orig.base;head=new complex[orig.size];
  93.    for(int i=0;i<size;i++)head[i]=orig.head[i];
  94.    }
  95.  
  96. Cmatrix operator*(Cmatrix&a, Cmatrix& b)
  97.    {
  98.    int i,j,k;complex zero;
  99.    if(a.colkt!=b.rowkt)
  100.       error(" cannot multiply matrices colkt,rowkt ",a.colkt,b.rowkt);
  101.    Cmatrix prod(a.rowkt,b.colkt,a.base);
  102.    for(i=0;i<a.rowkt;i++)for(j=0;j<b.colkt;j++)
  103.       {prod.m[i][j]=zero;
  104.       for(k=0;k<a.colkt;k++)prod.m[i][j] += a.m[i][k]*b.m[k][j];
  105.       }
  106.    return prod;
  107.    }   
  108.  
  109.  
  110. Cvector operator*(Cmatrix& a, Cvector& b)
  111.    {
  112.    int k;complex zero;
  113.    if(a.colkt!=b.size)error(" cannot mult. matrix,vector ",a.colkt,b.size);
  114.    Cvector prod(b.size,b.base,zero);
  115.    for(int i=0;i<a.rowkt;i++)
  116.       {prod.head[i]=zero;
  117.       for(k=0;k<a.colkt;k++)prod.head[i] += (a.m[i][k])*(b.head[k]);
  118.       }
  119.    return prod;
  120.    }   
  121.  
  122. complex Cvector::operator*(Cvector& rvalue) //dot product
  123.    {int i; complex *element,*relement; complex sum(0.,0.);// sum 
  124.    if(rvalue.size!= size)
  125.       error(" dot product unequal length vectors ",size,rvalue.size);
  126.    for(i=0,element=head,relement=rvalue.head;i< size;
  127.       i++,element++,relement++)
  128.       sum += *element * *relement;
  129.    return sum;
  130.    }
  131.  
  132. void Cvector::operator=(  Cvector& rhs)
  133.    {
  134.    if( this== &rhs)return;
  135.    if(size!=rhs.size)
  136.       {
  137.       delete head;
  138.       head= new complex[rhs.size];
  139.       }
  140.    for(int j=0;j<size;j++){head[j]=rhs.head[j];}
  141.    }
  142.  
  143. main() 
  144. {
  145. complex a,b(1.,0.),c(0.,1.); double dot; int i,j;
  146. printf(" size of complex=%d\n",sizeof(complex));
  147. a=b+c;
  148. a.print("summand=","\n");
  149. printf(" magnitude=%e\n", a.abs());
  150. b=a/c;
  151. b.print("quotient=","\n");
  152. a=b.cexp();
  153. a.print(" exponential=","\n");
  154. a= complex(1.,0.);b=complex(3.,2.);
  155. printf(" size of Cvector=%d\n",sizeof(Cvector));
  156. //
  157.    {Cvector A(4,0,a);Cvector B(4,0,b);Cvector D(4,0,b);
  158. for(j=0;j<B.size;j++) B.element(j).print(" B=","\n");   
  159.    a= A*B;
  160.    B=A;
  161. for(j=0;j<B.size;j++) B.element(j).print(" B=","\n");   
  162.    a.print(" dot product=","\n");
  163.       {Cmatrix m(4,4,0),n(4,4,0),o(4,4,0);
  164.       c=complex(0.,1.);
  165.       for (i=0;i<4;i++)
  166.          for(j=0;j<4;j++){
  167.             m.m[i][j].print(" matrix element=","\n");}
  168.       for (i=0;i<4;i++)
  169.          for(j=0;j<4;j++){
  170.             n.setelement( i,j, complex((double)(i+j),0.));
  171.             n.m[i][j].print(" matrix element=","\n");}
  172.       D = (m * A);
  173. for(j=0;j<B.size;j++) B.element(j).print(" B=","\n");   
  174.       B = (m * A);
  175. for(j=0;j<B.size;j++) B.element(j).print(" B=","\n");   
  176.       o= m * n;
  177.  
  178.       for (i=0;i<4;i++)
  179.          for(j=0;j<4;j++){
  180.             o.m[i][j].print(" matrix element=","\n");}
  181.       printf(" end of inner block\n");
  182.       }
  183. printf(" end of block\n");
  184.    }
  185. printf(" end of program\n");
  186. }
  187.  
  188.  
  189.