0000/*********************************************
0001
0002 Miscellaneous Functions
0003
0004********************************************/
0005#include <stdlib.h>
0006#include <math.h>
0007#include "gcs2.h"
0008
0009/********************************************
0010
0011 matric multiplication c=a*b
0012                   (l,m)=(l,n)*(n,m)
0013********************************************/
0014int matmul(darray *c,darray *a,darray *b)
0015{
0016 int i,j,k,l,m,n;
0017
0018 l=a->size[0];
0019 n=a->size[1];
0020 m=b->size[1];
0021 for(i=0;i<l;i++)
0022  for(j=0;j<m;j++){
0023   c->dat[i*m+j]=0.;
0024   for(k=0;k<n;k++)
0025    c->dat[i*m+j]+=a->dat[i*n+k]*b->dat[k*m+j];
0026  }
0027}
0028
0029/*******************************************
0030
0031 matrix inversion
0032
0033*******************************************/
0034int matinv(darray *a)
0035{
0036 int i,j,k,l,m,n;
0037 int irow,icolumn;
0038 double determ,amax,swap,t;
0039 int *ipivot;
0040 int *index1;
0041 int *index2;
0042 double *pivot;
0043
0044 n=a->size[0];
0045 ipivot=(int *)malloc(sizeof(int)*n);
0046 index1=(int *)malloc(sizeof(int)*n);
0047 index2=(int *)malloc(sizeof(int)*n);
0048 pivot=(double *)malloc(sizeof(double)*n);
0049 determ=1.;
0050 for(i=0;i<n;i++) ipivot[i]=0;
0051 for(i=0;i<n;i++){
0052/* search for pivot element */
0053  amax=0.;
0054  for(j=0;j<n;j++){
0055   if(ipivot[j]!=1){
0056    for(k=0;k<n;k++){
0057     if(ipivot[k]!=1){
0058      if(fabs(amax)<fabs(a->dat[j*n+k])){
0059       irow=j;
0060       icolumn=k;
0061       amax=a->dat[j*n+k];
0062      }
0063     }
0064    } /* for(k=0;k<n;k++) */
0065   }
0066  } /* for(j=0;j<n;j++) */
0067  ipivot[icolumn]=1;
0068/* interchange rows to put pivot element on diagonal */
0069  if((irow+icolumn)%2==1) determ=-determ;
0070  if(irow!=icolumn){
0071   for(l=0;l<n;l++){
0072    swap=a->dat[irow*n+l];
0073    a->dat[irow*n+l]=a->dat[icolumn*n+l];
0074    a->dat[icolumn*n+l]=swap;
0075   }
0076  }
0077  index1[i]=irow;
0078  index2[i]=icolumn;
0079  pivot[i]=a->dat[icolumn*n+icolumn];
0080  determ*=pivot[i];
0081/* divide pivot row by pivot element */
0082  a->dat[icolumn*n+icolumn]=1.;
0083  for(l=0;l<n;l++) a->dat[icolumn*n+l]/=pivot[i];
0084/* reduce non-pivot rows */
0085  for(j=0;j<n;j++){
0086   if(j!=icolumn){
0087    t=a->dat[j*n+icolumn];
0088    a->dat[j*n+icolumn]=0.;
0089    for(k=0;k<n;k++) a->dat[j*n+k]-=a->dat[icolumn*n+k]*t;
0090   }
0091  }
0092 } /* for(i=0;i<n;i++) */
0093/* interchange columns */
0094 for(i=n-1;i>=0;i--){
0095  if(index1[i]!=index2[i]){
0096   irow=index1[i];
0097   icolumn=index2[i];
0098   for(k=0;k<n;k++){
0099    swap=a->dat[k*n+irow];
0100    a->dat[k*n+irow]=a->dat[k*n+icolumn];
0101    a->dat[k*n+icolumn]=swap;
0102   }
0103  }
0104 }
0105 free(ipivot);
0106 free(index1);
0107 free(index2);
0108 free(pivot);
0109 return;
0110}
0111