
/****Copyright (c) 2001  Alexandros V. Gerbessiotis. 
 **** All rights reserved.                               
 ****Permission to use, copy, modify, and distribute this software,
 ****and to incorporate it, in whole or in part, into 
 ****other software,is hereby granted without fee, provided that
 ****  (1) the above copyright notice and this permission 
 ****      notice appear in all copies of the source code, and 
 ****      the above copyright notice appear in clearly visible 
 ****      form on all supporting documentation and distribution 
 ****      media, and
 ****     (2) any redistribution of the software, in original
 ****         or modified form, be without fee and subject to
 ****         these same conditions.
 ****   No guarantee, whatsoever, is offered that the code works as
 ****   advertised or is absent of any, even damaging,
 ****   side-effects. Use this code at your own personal risk.
 ****/
#include	<math.h>

void transpose(double *a, double *b, int n)
{
 register int i,j,k,l;
 register double temp;

 for(i=0;i<n;i++) {
   k=n*i;l=k;
   for(j=i;j<n;j++) {
      temp= a[k  +j];	
      a[k  +j]=b[l  +i];
      b[l  +i] = temp;
      l+=n;
   }
 }
}

void fast(double *a,int n) {
 register int i,j,m;
 register double temp;
 m=n*n-2;
 temp=0.0;
 for(j=0;j<n;j++) {
   for(i=0;i<m;i++) {
      temp += a[i] * a[i+1];
   }
 a[0] = temp;
 }
}
void mult3(double *a,double *b, double *c, int n) {
register int i,j,k,l,m;
register double temp;

 for(j=0,m=0;j<n;++j,m+=n){
   for(i=0,l=0;i<n;++i,l+=n) {
     temp=0.0;
     for(k=0;k<n;++k)
       temp += a[l+k] * b[m+k] ;
     c[m+i] = temp;
   }
 }
}


void mult(double *a,double *b, double *c, int n) {
register int i,j,k,l,m;
register double temp;

 for(i=0,l=0;i<n;++i,l+=n)
   for(j=0,m=0;j<n;++j,m+=n) {
     temp=0.0;
     for(k=0;k<n;++k)
       temp += a[l+k] * b[m+k] ;
     c[m+i] = temp;
   }
}

void mult2(double *a,double *b, double *c, int n) {
 register int i,j,k,l,m;
 register double temp;

 for(j=0;j<n;++j){
   m=j*n;
   for(i=0;i<n;++i) {
     l=i*n;
     temp=0.0;
     for(k=0;k<n;++k)
       temp += a[l+k] * b[m+k] ;
     c[m+i] = temp;
   }
 }
}
