utility.c

矩阵运算的有关部分代码 不全 为了注册才上载的 以后有好东东了在给各位大虾分享

// -*- C -*-
// Werner Stuetzle, 1-16-03
//-------------------------

#include<assert.h>
#include<stdio.h>
#include<stdlib.h>
//#include "nrutil.h"
//#include "least-squares.h"
#include<math.h>


void matrix_times_matrix (double **a, double **b, double **result,
		      int nrow_a, int ncol_a, int ncol_b) {
  int i, j, k;
  double sum;
  for (i = 0; i < nrow_a; i++) {
    for (j = 0; j < ncol_b; j++) {
      sum = 0.0;
      for (k =0; k < ncol_a; k++) {
	     sum = sum + a[i][k] * b[k][j];
      }
      result[i][j] = sum;
    }
  }
}

/*----------------------------------------------------------------------*/

void matrix_times_vector (double **a, double *x, double *result,
			  int nrow, int ncol) {
  int i, j;
  double sum;
  for (i = 0; i < nrow; i++) {
    sum = 0.0;
    for (j = 0; j < ncol; j++) {
      sum = sum + a[i][j] * x[j];
    }
    result[i] = sum;
  }
}

/*----------------------------------------------------------------------*/

void print_vector (double *a, int length, FILE *stream ) {
  int i;
  for (i = 0; i < length; i++) {
    fprintf(stream, "%f   ", a[i]);
  }
  fprintf(stream, "\n");
}

/*----------------------------------------------------------------------*/

void print_matrix (double **a, int nrow, int ncol, FILE* stream) {
  int i;
  for (i = 0; i < nrow; i++) {
    print_vector (a[i], ncol, stream);
  }
}

/*---------------------------------------------------------------------*/

void copy_matrix (double **from, double **to, int nrow, int ncol) {
  int i, j;
  for (i = 0; i < nrow; i++) {
    for (j = 0; j < ncol; j++) {
      to[i][j] = from[i][j];
    }
  }
}

/*----------------------------------------------------------------------*/

void transpose_matrix (double **a, double **a_transpose, int nrow, int ncol) {
  int i, j;
  for (i = 0; i < nrow; i++) {
    for (j = 0; j < ncol; j++) {
      a_transpose[j][i] = a[i][j];
    }
  }
}

static void SWAP(double *a, double *b)
{
        double temp;
        temp = *a;
        *a = *b;
        *b = temp;
}

static void InvMatrix(double **a, int n)
{
        int *indxc,*indxr,*ipiv;
        int i,icol,irow,j,k,l,ll;
        double big,dum,pivinv,temp;

        indxc=ivector(0,n-1);
        indxr=ivector(0,n-1);
        ipiv=ivector(0,n-1);

        for (j=0;j<=n-1;j++) 
                ipiv[j]=0;

        for (i=0;i<=n-1;i++)
        {
                big=0.0;
                for (j=0;j<=n-1;j++)
                        if (ipiv[j] != 1)
                                for (k=0;k<=n-1;k++)
                                {
                                        if (ipiv[k] == 0)
                                        {
                                                if (fabs(a[j][k]) >= big)
                                                {
                                                        big=fabs(a[j][k]);
                                                        irow=j;
                                                        icol=k;
                                                }
                                        } 
                                        else if (ipiv[k] > 1)
                                                nrerror("Inverse Matrix: Singular Matrix-1");
                                }
                ++(ipiv[icol]);
                if (irow != icol) 
                {
                        for (l=0;l<=n-1;l++)
                                SWAP(&a[irow][l],&a[icol][l]);
                }
                indxr[i]=irow;
                indxc[i]=icol;
                if (a[icol][icol] == 0.0) 
                        nrerror("Inverse Matrix: Singular Matrix-2");
                pivinv=1.0/a[icol][icol];
                a[icol][icol]=1.0;
                for (l=0;l<=n-1;l++)
                        a[icol][l] *= pivinv;
                
                for (ll=0;ll<=n-1;ll++)
                        if (ll != icol) 
                        {
                                dum=a[ll][icol];
                                a[ll][icol]=0.0;
                                for (l=0;l<=n-1;l++)
                                        a[ll][l] -= a[icol][l]*dum;
                        }
        }

        for (l=n-1;l>=0;l--) 
        {
                if (indxr[l] != indxc[l])
                        for (k=0;k<=n-1;k++)
                                SWAP(&a[k][indxr[l]],&a[k][indxc[l]]);
        }

        free_ivector(ipiv,0,n-1);
        free_ivector(indxr,0,n-1);
        free_ivector(indxc,0,n-1);
}