matinvert.cpp

模糊聚類分析源碼。包含教學文件

/*
  
  Context       : Matrix and Vector Operation

  Author        : Frank Hoeppner, see also AUTHORS file 

  Description   : implementation of function module matinvert
                  
  History       :
    matinvert.nw:
     990218 fh: first version

  Comment       : 
    This file was generated automatically. DO NOT EDIT.

  Copyright     : Copyright (C) 1999-2000 Frank Hoeppner

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#ifndef matinvert_SOURCE
#define matinvert_SOURCE

/* configuration include */
#ifdef HAVE_CONFIG_H
/*
//FILETREE_IFDEF HAVE_CONFIG_H
*/
#include "config.h"
/*
//FILETREE_ENDIF
*/
#endif



/* necessary includes */
#include "matinvert.hpp"
#include "trace.hpp" // Log/Trace

/* private typedefs */


/* private functions */


/* data */
const char MATINVERT_ID[] = "MatInvert";

/* implementation */
template <class M>
void gauss_jordan(M& A)
  {
  FUNCTION_LOG(MATINVERT_ID,"gauss_jordan");
  
  int hi,i,j,k,r,pivot[A.cols()],n = A.cols();
  typename M::value_type max,hr,hv[A.cols()];
  
  invariant( n==A.rows(), "inv(AbsMatrix) (dimension mismatch)");
  //invariant( n<10 , "inv(AbsMatrix) (fixed range exceeded)");
  
  // Gau\3-Jordan-Verfahren, Stoer S. 161
  for (j=0;j<n;++j) 
    {
    pivot[j]=j;
    }
    
  for (j=0;j<n;++j)
    {
    // Pivotsuche (bestimme Pivotelement fuer j.-te Zeile
    max = fabs(A(j,j)); 
    r = j;
    trace(MATINVERT_ID,"search for pivot element",make_pair(j,A));
    for (i=j+1;i<n;++i)
      {
      hr = fabs(A(i,j));
      if (hr>max) 
        { 
        max=hr; 
        r=i; 
        }
      }
    trace(MATINVERT_ID,"select row",r);
    invariant(max!=0.0,"gauss_jordan (singularity)");

    // Zeilentausch
    if (r>j)
      {
      for (k=0;k<n;k++) 
        { 
        hr=A(j,k); 
        A(j,k)=A(r,k); 
        A(r,k)=hr;     
        }

      hi=pivot[j]; 
      pivot[j]=pivot[r]; 
      pivot[r]=hi;
      }


    // Transformation
    hr=1.0/A(j,j);
    for (i=0;i<n;i++) A(i,j)*=hr;
    A(j,j)=hr;
    for (k=0;k<n;k++) if (k!=j)
      {
      for (i=0;i<n;i++) if (i!=j)
        A(i,k)-=A(i,j)*A(j,k);
      A(j,k)*=-hr;
      }
    }

  // Spaltentausch
  for (i=0;i<n;i++)
    {
    for (k=0;k<n;k++) hv[pivot[k]]=A(i,k);
    for (k=0;k<n;k++) A(i,k)=hv[k];
    }
  }

/* template instantiation */
#include "Matrix.hpp"
#ifdef _SGI_SOURCE
#  pragma instantiate void gauss_jordan(Matrix<4,4>&)
#  pragma instantiate void gauss_jordan(Matrix<3,3>&)
#endif
#ifdef __GNUC__
//   template void gauss_jordan(Matrix<4,4>&);
//   template void gauss_jordan(Matrix<3,3>&);
#endif


#endif /* matinvert_SOURCE */