The C Users' Group Library 1994 August/wc-cdrom-cusersgrouplibrary-1994-08.iso/listings/v_09_01/9n01128a

Text File | 1990-03-26 | 5KB | 189 lines

#include <stdlib.h>
#include "cmatrix.hpp"

error(char* msg,int index,int s)
   {
   printf("%s%d %d\n",msg,index,s);
   exit(1);
   }

void Cvector::check( int index)
   {int loc;
   loc=index-(*this).base;
   if (loc<0)
      error(" Cvector error: index-base<0",index,(*this).base);
   if (loc>= (*this).size)
      error(" Cvector error: index too large",index,(*this).size);
   }

void Cmatrix::check( int i,int j)
   {int loc;
   loc=i-(*this).base;
   if (loc<0)
      error(" Cmatrix error: row index-base<0",i,(*this).base);
   if (loc>= (*this).rowkt)
      error(" Cmatrix error: row index too large",i,(*this).rowkt);

   loc=j-(*this).base;
   if (loc<0)
      error(" Cmatrix error: column index-base<0",j,(*this).base);
   if (loc>= (*this).colkt)
      error(" Cmatrix error: column index too large",j,(*this).colkt);

   }


void Cmatrix::init(int row, int col, int b)
   {
   if(row<1||col<1)error(" matrix needs at least one row/col ",row,col);
   rowkt=row;colkt=col;base=b;
   m= new complex *[rowkt];
   for(int j=0;j<rowkt;j++)m[j]= new complex[colkt];
   complex zero;complex one(1.,0.);
   //initialize to identity matrix.
   for(int i=0;i<rowkt;i++){for(j=0;j<colkt;j++)m[i][j]=zero;
      m[i][i]= one;// omit this line for init. to zero matrix
      }
   }

Cmatrix::Cmatrix(Cmatrix& a)
   { init(a.rowkt,a.colkt,a.base);
    for(int i=0;i<a.rowkt;i++)for(int j=0;j<a.colkt;j++)m[i][j]=a.m[i][j];
   }

Cmatrix::~Cmatrix()
   {for(int i=0;i<rowkt;i++)delete m[i];
   delete m;
   }


complex& Cvector::element(int i)
   {
   check(i);return head[i-base];
   }

void Cvector::setelement( int i, complex& value)
   {
   check(i);
   head[i]=value;
   }
complex& Cmatrix::element(int i,int j)
   {
   check(i,j);return m[i-base][j-base];
   }

void Cmatrix::setelement( int i, int j, complex value)
   {
   check(i,j);m[i][j]=value;
   }


void Cmatrix::operator=(Cmatrix& a)
   {
   if(this==&a)return;
   for(int i=0;i<rowkt;i++)delete m[i];
   delete m;
   init(a.rowkt,a.colkt,a.base);
   for(i=0;i<a.rowkt;i++)for(int j=0;j<a.colkt;j++)m[i][j]=a.m[i][j];
   }

 Cvector::Cvector( Cvector& orig) //copy
   {
   size=orig.size;base=orig.base;head=new complex[orig.size];
   for(int i=0;i<size;i++)head[i]=orig.head[i];
   }

Cmatrix operator*(Cmatrix&a, Cmatrix& b)
   {
   int i,j,k;complex zero;
   if(a.colkt!=b.rowkt)
      error(" cannot multiply matrices colkt,rowkt ",a.colkt,b.rowkt);
   Cmatrix prod(a.rowkt,b.colkt,a.base);
   for(i=0;i<a.rowkt;i++)for(j=0;j<b.colkt;j++)
      {prod.m[i][j]=zero;
      for(k=0;k<a.colkt;k++)prod.m[i][j] += a.m[i][k]*b.m[k][j];
      }
   return prod;
   }


Cvector operator*(Cmatrix& a, Cvector& b)
   {
   int k;complex zero;
   if(a.colkt!=b.size)error(" cannot mult. matrix,vector ",a.colkt,b.size);
   Cvector prod(b.size,b.base,zero);
   for(int i=0;i<a.rowkt;i++)
      {prod.head[i]=zero;
      for(k=0;k<a.colkt;k++)prod.head[i] += (a.m[i][k])*(b.head[k]);
      }
   return prod;
   }

complex Cvector::operator*(Cvector& rvalue) //dot product
   {int i; complex *element,*relement; complex sum(0.,0.);// sum
   if(rvalue.size!= size)
      error(" dot product unequal length vectors ",size,rvalue.size);
   for(i=0,element=head,relement=rvalue.head;i< size;
      i++,element++,relement++)
      sum += *element * *relement;
   return sum;
   }

void Cvector::operator=(  Cvector& rhs)
   {
   if( this== &rhs)return;
   if(size!=rhs.size)
      {
      delete head;
      head= new complex[rhs.size];
      }
   for(int j=0;j<size;j++){head[j]=rhs.head[j];}
   }

main()
{
complex a,b(1.,0.),c(0.,1.); double dot; int i,j;
printf(" size of complex=%d\n",sizeof(complex));
a=b+c;
a.print("summand=","\n");
printf(" magnitude=%e\n", a.abs());
b=a/c;
b.print("quotient=","\n");
a=b.cexp();
a.print(" exponential=","\n");
a= complex(1.,0.);b=complex(3.,2.);
printf(" size of Cvector=%d\n",sizeof(Cvector));
//
   {Cvector A(4,0,a);Cvector B(4,0,b);Cvector D(4,0,b);
for(j=0;j<B.size;j++) B.element(j).print(" B=","\n");
   a= A*B;
   B=A;
for(j=0;j<B.size;j++) B.element(j).print(" B=","\n");
   a.print(" dot product=","\n");
      {Cmatrix m(4,4,0),n(4,4,0),o(4,4,0);
      c=complex(0.,1.);
      for (i=0;i<4;i++)
         for(j=0;j<4;j++){
            m.m[i][j].print(" matrix element=","\n");}
      for (i=0;i<4;i++)
         for(j=0;j<4;j++){
            n.setelement( i,j, complex((double)(i+j),0.));
            n.m[i][j].print(" matrix element=","\n");}
      D = (m * A);
for(j=0;j<B.size;j++) B.element(j).print(" B=","\n");
      B = (m * A);
for(j=0;j<B.size;j++) B.element(j).print(" B=","\n");
      o= m * n;

      for (i=0;i<4;i++)
         for(j=0;j<4;j++){
            o.m[i][j].print(" matrix element=","\n");}
      printf(" end of inner block\n");
      }
printf(" end of block\n");
   }
printf(" end of program\n");
}