yimop.h

[[ Complex Matrices : Language c]]

/* xbfract.h                    free ware                  [e-mail]  */


/* --------------------------------- FUNCTION ------------------------------  */
/* Do    : mA + mB = mAplsB                                                   */
/*                                                                            */
/* -------------------------------------------------------------------------- */
void add_mzF(
PmzF mA,
PmzF mB,
PmzF mAplsB)
{
     int i;
     int j;
complexF za;
complexF zb;
complexF zaplsb;

  if (mA->rows != mB->rows || mA->cols != mB->cols)
    {
     printf("\n addm() error - matrices different sizes");
     printf("\n Press return to continue");
     getchar();
     exit(1);
    }
  if (mA->rows!=mAplsB->rows||mA->cols!=mAplsB->cols)
    {
     printf("\n addm() error - matrices different sizes");
     printf("\n Press return to continue");
     getchar();
     exit(1);
    }
  for (   i = 0; i < mA->rows ; i++)
    for ( j = 0; j < mA->cols ; j += C4)
         {
          za.r.n = *(mA->pb+i *mA->cols+j  );
          za.r.d = *(mA->pb+i *mA->cols+j+1);
          za.i.n = *(mA->pb+i *mA->cols+j+2);
          za.i.d = *(mA->pb+i *mA->cols+j+3);

          zb.r.n = *(mB->pb+i *mB->cols+j  );
          zb.r.d = *(mB->pb+i *mB->cols+j+1);
          zb.i.n = *(mB->pb+i *mB->cols+j+2);
          zb.i.d = *(mB->pb+i *mB->cols+j+3);

          zaplsb = add_zF(za, zb);

          *(mAplsB->pb+i *mAplsB->cols+j  )   = zaplsb.r.n;
          *(mAplsB->pb+i *mAplsB->cols+j+1)   = zaplsb.r.d;
          *(mAplsB->pb+i *mAplsB->cols+j+2)   = zaplsb.i.n;
          *(mAplsB->pb+i *mAplsB->cols+j+3)   = zaplsb.i.d;
         }
}


/* --------------------------------- FUNCTION ------------------------------  */
/* Do    : mA + mB = mAmnsB                                                   */
/*                                                                            */
/* -------------------------------------------------------------------------- */
void sub_mzF(
PmzF mA,
PmzF mB,
PmzF mAmnsB)
{
     int i;
     int j;
complexF za;
complexF zb;
complexF zamnsb;

  if (mA->rows != mB->rows || mA->cols != mB->cols)
    {
     printf("\n addm() error - matrices different sizes");
     printf("\n Press return to continue");
     getchar();
     exit(1);
    }
  if (mA->rows!=mAmnsB->rows||mA->cols!=mAmnsB->cols)
    {
     printf("\n addm() error - matrices different sizes");
     printf("\n Press return to continue");
     getchar();
     exit(1);
    }
  for (   i = 0; i < mA->rows ; i++)
    for ( j = 0; j < mA->cols ; j += C4)
         {
          za.r.n = *(mA->pb+i *mA->cols+j  );
          za.r.d = *(mA->pb+i *mA->cols+j+1);
          za.i.n = *(mA->pb+i *mA->cols+j+2);
          za.i.d = *(mA->pb+i *mA->cols+j+3);

          zb.r.n = *(mB->pb+i *mB->cols+j  );
          zb.r.d = *(mB->pb+i *mB->cols+j+1);
          zb.i.n = *(mB->pb+i *mB->cols+j+2);
          zb.i.d = *(mB->pb+i *mB->cols+j+3);

          zamnsb = sub_zF(za, zb);

          *(mAmnsB->pb+i *mAmnsB->cols+j  )   = zamnsb.r.n;
          *(mAmnsB->pb+i *mAmnsB->cols+j+1)   = zamnsb.r.d;
          *(mAmnsB->pb+i *mAmnsB->cols+j+2)   = zamnsb.i.n;
          *(mAmnsB->pb+i *mAmnsB->cols+j+3)   = zamnsb.i.d;
         }
}


/* --------------------------------- FUNCTION ------------------------------  */
/* Do    : mA * mB = mAB                                                      */
/*                                                                            */
/* -------------------------------------------------------------------------- */
void mul_mzF(
PmzF mA,
PmzF mB,
PmzF mAB)
{
     int i;
     int j;
     int k;
complexF za;
complexF zb;
complexF zab;
complexF zabplsab;

    for( k = 0; k < mA->rows ; k++)
     for( j = 0; j < mB->cols ; j += C4)
     {
      zabplsab = initD_zF(0,1,  0,1);
      for( i = 0; i < mA->cols ; i += C4)
         {
          za.r.n = *(mA->pb+k *mA->cols+i  );
          za.r.d = *(mA->pb+k *mA->cols+i+1);
          za.i.n = *(mA->pb+k *mA->cols+i+2);
          za.i.d = *(mA->pb+k *mA->cols+i+3);

          zb.r.n = *(mB->pb+(i/C4) *mB->cols+j  );
          zb.r.d = *(mB->pb+(i/C4) *mB->cols+j+1);
          zb.i.n = *(mB->pb+(i/C4) *mB->cols+j+2);
          zb.i.d = *(mB->pb+(i/C4) *mB->cols+j+3);

          zab = mul_zF(za,zb);

          zabplsab = add_zF(zabplsab,zab);
         }
         *(mAB->pb+k *mAB->cols+j  ) = zabplsab.r.n;
         *(mAB->pb+k *mAB->cols+j+1) = zabplsab.r.d;
         *(mAB->pb+k *mAB->cols+j+2) = zabplsab.i.n;
         *(mAB->pb+k *mAB->cols+j+3) = zabplsab.i.d;
     }
}


/* --------------------------------- FUNCTION ------------------------------  */
/* Do    : mA**n = mAPn                                                       */
/*                                                                            */
/* -------------------------------------------------------------------------- */
void pow_mzF(
PmzF mA,
int     Pn,
PmzF mAPn
)
{
 int i;

double pbT[MXR][MXC*C4];mzF mT ={MXR,MXC*C4,&pbT[0][0]};

            mT.rows = mA->rows;
            mT.cols = mA->rows*C4;

 copy_mzF(mA,mAPn);

 if(!Pn) mid_mzF(mAPn);
 else    for(i = Pn-1 ; i ; i--)
            {
             mul_mzF(mAPn,mA,&mT);
             copy_mzF(&mT,mAPn);
            }
}


/* --------------------------------- FUNCTION ------------------------------  */
/* Do    :                                                                    */
/*                                                                            */
/* -------------------------------------------------------------------------- */
void Fmul_mzF(
fraction f,
PmzF  mA,
PmzF  mfA)
{
     int i;
     int j;
complexF z;

  for (   i = 0; i < mA->rows ; i++)
    for ( j = 0; j < mA->cols ; j += C4)
         {
          z.r.n = *(mA->pb+i *mA->cols+j  );
          z.r.d = *(mA->pb+i *mA->cols+j+1);
          z.i.n = *(mA->pb+i *mA->cols+j+2);
          z.i.d = *(mA->pb+i *mA->cols+j+3);

          z = Fmul_zF(f,z);

          *(mfA->pb+i *mfA->cols+j  ) = z.r.n;
          *(mfA->pb+i *mfA->cols+j+1) = z.r.d;
          *(mfA->pb+i *mfA->cols+j+2) = z.i.n;
          *(mfA->pb+i *mfA->cols+j+3) = z.i.d;
         }
}


/* --------------------------------- FUNCTION ------------------------------  */
/* Do    :                                                                    */
/*                                                                            */
/* -------------------------------------------------------------------------- */
void Zmul_mzF(
complexF za,
PmzF  mA,
PmzF  mfA)
{
     int i;
     int j;
complexF z;

  for (   i = 0; i < mA->rows ; i++)
    for ( j = 0; j < mA->cols ; j += C4)
         {
          z.r.n = *(mA->pb+i *mA->cols+j  );
          z.r.d = *(mA->pb+i *mA->cols+j+1);
          z.i.n = *(mA->pb+i *mA->cols+j+2);
          z.i.d = *(mA->pb+i *mA->cols+j+3);

          z = mul_zF(za,z);

          *(mfA->pb+i *mfA->cols+j  ) = z.r.n;
          *(mfA->pb+i *mfA->cols+j+1) = z.r.d;
          *(mfA->pb+i *mfA->cols+j+2) = z.i.n;
          *(mfA->pb+i *mfA->cols+j+3) = z.i.d;
         }
}

/* --------------------------------- FUNCTION ------------------------------  */
/* Do    :                                                                    */
/*                                                                            */
/* -------------------------------------------------------------------------- */
void cnj_mzF(
PmzF mA,
PmzF mConjA)
{
     int i;
     int j;
complexF z;

  for (   i = 0; i < mA->rows ; i++)
    for ( j = 0; j < mA->cols ; j += C4)
         {
          z.r.n = *(mA->pb+i *mA->cols+j  );
          z.r.d = *(mA->pb+i *mA->cols+j+1);
          z.i.n = *(mA->pb+i *mA->cols+j+2);
          z.i.d = *(mA->pb+i *mA->cols+j+3);

          z = cnj_zF(z);

                  *(mConjA->pb+i *mConjA->cols+j  ) = z.r.n;
                  *(mConjA->pb+i *mConjA->cols+j+1) = z.r.d;
                  *(mConjA->pb+i *mConjA->cols+j+2) = z.i.n;
                  *(mConjA->pb+i *mConjA->cols+j+3) = z.i.d;
         }
}


/* --------------------------------- FUNCTION ------------------------------  */
/* Do    :                                                                    */
/*                                                                            */
/* -------------------------------------------------------------------------- */
void cnj_transpose_mzF(
PmzF mA,
PmzF mTrpsA)
{
     int i;
     int j;
complexF z;

  if (mTrpsA->rows * C4 != mA->cols)
    {
     printf("\n transpose() error - dest matrix rows must = source matrix cols");
     printf("\n Press return to continue");
     getchar();
     exit(1);
    }
  if (mA->rows * C4 != mTrpsA->cols)
    {
     printf("\n transpose() error - source matrix rows must = dest matrix cols");
     printf("\n Press return to continue");
     getchar();
     exit(1);
    }

  for (   i = 0; i < mA->rows ; i++)
    for ( j = 0; j < mA->cols ; j += C4)
         {
          z.r.n = *(mA->pb+i *mA->cols+j  );
          z.r.d = *(mA->pb+i *mA->cols+j+1);
          z.i.n = *(mA->pb+i *mA->cols+j+2);
          z.i.d = *(mA->pb+i *mA->cols+j+3);

          z = cnj_zF(z);

          *(mTrpsA->pb+j/C4 *mTrpsA->cols+i*C4  )   = z.r.n;
          *(mTrpsA->pb+j/C4 *mTrpsA->cols+i*C4+1)   = z.r.d;
          *(mTrpsA->pb+j/C4 *mTrpsA->cols+i*C4+2)   = z.i.n;
          *(mTrpsA->pb+j/C4 *mTrpsA->cols+i*C4+3)   = z.i.d;
         }
}


/* --------------------------------- FUNCTION ------------------------------  */
/* Do    :                                                                    */
/*                                                                            */
/* -------------------------------------------------------------------------- */
void transpose_mzF(
PmzF mA,
PmzF mTrpsA)
{
     int i;
     int j;
complexF z;

  if (mTrpsA->rows * C4 != mA->cols)
    {
     printf("\n transpose() error - dest matrix rows must = source matrix cols");
     printf("\n Press return to continue");
     getchar();
     exit(1);
    }
  if (mA->rows * C4 != mTrpsA->cols)
    {
     printf("\n transpose() error - source matrix rows must = dest matrix cols");
     printf("\n Press return to continue");
     getchar();
     exit(1);
    }

  for (   i = 0; i < mA->rows ; i++)
    for ( j = 0; j < mA->cols ; j += C4)
         {
          z.r.n = *(mA->pb+i *mA->cols+j  );
          z.r.d = *(mA->pb+i *mA->cols+j+1);
          z.i.n = *(mA->pb+i *mA->cols+j+2);
          z.i.d = *(mA->pb+i *mA->cols+j+3);

          *(mTrpsA->pb+j/C4 *mTrpsA->cols+i*C4  )   = z.r.n;
          *(mTrpsA->pb+j/C4 *mTrpsA->cols+i*C4+1)   = z.r.d;
          *(mTrpsA->pb+j/C4 *mTrpsA->cols+i*C4+2)   = z.i.n;
          *(mTrpsA->pb+j/C4 *mTrpsA->cols+i*C4+3)   = z.i.d;
         }
}

/* --------------------------------- FUNCTION ------------------------------  */
/* Do    :                                                                    */
/*                                                                            */
/* -------------------------------------------------------------------------- */
complexF trace_mzF(
PmzF m)
{
     int i;
complexF Trace;
complexF zT;

  if ((m->rows * C4) !=  m->cols)
    {
     printf("\n trF() error - Square mF, please");
     printf("\n Press return to continue");
     getchar();
     exit(1);
    }

    Trace = initD_zF(0,1,0,1);

  for(   i = 0; i < m->rows ; i++)
     {
      zT.r.n = (*(m->pb+i *m->cols+i*C4   ));
      zT.r.d = (*(m->pb+i *m->cols+i*C4+1));
      zT.i.n = (*(m->pb+i *m->cols+i*C4+2));
      zT.i.d = (*(m->pb+i *m->cols+i*C4+3));

             zT.r = mini_F(  zT.r);
             zT.i = mini_F(  zT.i);

           Trace  = add_zF(Trace,zT);

         Trace.r  = mini_F(Trace.r);
         Trace.i  = mini_F(Trace.i);
     }

  return(Trace);
}