matrix.java

化学图形处理软件

/* $RCSfile$
 * $Author: egonw $
 * $Date: 2007-01-04 18:46:10 +0100 (Thu, 04 Jan 2007) $
 * $Revision: 7636 $
 *
 * Copyright (C) 2001-2007  The Chemistry Development Kit (CDK) project
 * 
 * Contact: [e-mail]
 * 
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1
 * of the License, or (at your option) any later version.
 * All we ask is that proper credit is given for our work, which includes
 * - but is not limited to - adding the above copyright notice to the beginning
 * of your source code files, and to any copyright notice that you may distribute
 * with programs based on this work.
 * 
 * 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 Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
 *
 */
package org.openscience.cdk.math;
 
import java.text.DecimalFormat;

/**
 * This class contains a matrix.
 *
 * @author  Stephan Michels <[e-mail]>
 * @cdk.created 2001-06-07
 */
public class Matrix
{
  // Attention! Variables are unprotected
  /** the content of this matrix **/
  public double[][] matrix;

  /** the number of rows of this matrix */
  public int rows;
  /** the number of columns of this matrix */
  public int columns;

  /**
   * Creates a new Matrix.
   */ 
  public Matrix(int rows, int columns)
  {
    this.rows = rows;
    this.columns = columns;
    matrix = new double[rows][columns];
  }

  /**
   * Creates a Matrix with content of an array.
   */
  public Matrix(double[][] array)
  {
    rows = array.length;
    int i,j;
    columns = array[0].length;
    for(i=1; i<rows; i++)
      columns = Math.min(columns,array[i].length);
    
    matrix = new double[rows][columns];
    for(i=0; i<rows; i++)
      for(j=0; j<columns; j++)
        matrix[i][j] = array[i][j];
  }

  /**
   * Returns the number of rows.
   */ 
  public int getRows()
  {
    return rows;
  }

  /**
   * Returns the number of columns.
   */
  public int getColumns()
  {
    return columns;
  }

  /**
   * Creates a Vector with the content of a row from this Matrix.
   */
  public Vector getVectorFromRow(int index)
  {
    double[] row = new double[columns];
    for(int i=0; i<columns; i++)
      row[i] = matrix[index][i];
    return new Vector(row);
  }

  /**
   * Creates a Vector with the content of a column from this Matrix.
   */
  public Vector getVectorFromColumn(int index)
  {
    double[] column = new double[rows];
    for(int i=0; i<rows; i++)
      column[i] = matrix[i][index];
    return new Vector(column);
  }

  /**
   * Creates a Vector with the content of the diagonal elements from this Matrix.
   */
  public Vector getVectorFromDiagonal()
  {
    int size = Math.min(rows, columns);
    Vector result = new Vector(size);
    for(int i=0; i<rows; i++)
      result.vector[i] = matrix[i][i];
    return result;
  }

  /**
   * Adds two matrices.
   */
  public Matrix add(Matrix b)
  {
    if ((b==null) ||
        (rows!=b.rows) || (columns!=b.columns))
      return null;
      
    int i, j;
    Matrix result = new Matrix(rows, columns);
    for(i=0; i<rows; i++)
      for(j=0; j<columns; j++)
        result.matrix[i][j] = matrix[i][j]+b.matrix[i][j];
    return result;
  }
  
  /**
   * Subtracts from two matrices.
   */
  public Matrix sub(Matrix b)
  {
    if ((b==null) ||
        (rows!=b.rows) || (columns!=b.columns))
      return null;
      
    int i, j;
    Matrix result = new Matrix(rows, columns);
    for(i=0; i<rows; i++) 
      for(j=0; j<columns; j++)
        result.matrix[i][j] = matrix[i][j]-b.matrix[i][j];
    return result;
  }

  /**
   * Multiplies this Matrix with another one.
   */
  public Matrix mul(Matrix b)
  {
    if ((b==null) ||
        (columns!=b.rows))
      return null;

    Matrix result = new Matrix(rows, b.columns);
    int i,j,k;
    double sum;
    for(i=0; i<rows; i++)
      for(k=0; k<b.columns; k++)
      {
        sum = 0;
        for(j=0; j<columns; j++)
          sum += matrix[i][j]*b.matrix[j][k];
        result.matrix[i][k] = sum;
      }

    return result;
  }
  
  /**
   *  Multiplies a Vector with this Matrix.
   */
  public Vector mul(Vector a)
  {
    if ((a==null) ||
        (columns!=a.size))
      return null;

    Vector result = new Vector(rows);
    int i,j;
    double sum;
    for(i=0; i<rows; i++)
    {
      sum = 0;
      for(j=0; j<columns; j++)
        sum += matrix[i][j]*a.vector[j];
      result.vector[i] = sum;
    }
    return result;
  }

  /**
   * Multiplies a scalar with this Matrix.
   */
  public Matrix mul(double a)
  {
    Matrix result = new Matrix(rows, columns);
    int i,j;
    for(i=0; i<rows; i++)
      for(j=0; j<columns; j++)
        result.matrix[i][j] = matrix[i][j]*a;

    return result;
  }

  /**
   * Copies a matrix.
   */
  public Matrix duplicate()
  {
    Matrix result = new Matrix(rows, columns);
    int i,j;
    for(i=0; i<rows; i++)
      for(j=0; j<columns; j++)
        result.matrix[i][j] = matrix[i][j];
    return result;
  }

  /**
   * Transposes a matrix.
   */
  public Matrix transpose()
  { 
    Matrix result = new Matrix(columns, rows);
    int i,j; 
    for(i=0; i<rows; i++) 
      for(j=0; j<columns; j++)
        result.matrix[j][i] = matrix[i][j];
    return result;
  }

  /**
   * Similar transformation
   * Ut * M * U
   */
  public Matrix similar(Matrix U)
  {
    Matrix result = new Matrix(U.columns, U.columns);
    double sum, innersum;
    for(int i=0; i<U.columns; i++)
      for(int j=0; j<U.columns; j++)
      {
        sum = 0d;
        for(int k=0; k<U.columns; k++)
        {
          innersum = 0d;
          for(int l=0; l<U.columns; l++)
            innersum += matrix[k][l]*U.matrix[l][j];
          sum += U.matrix[k][i]*innersum;
        }
        result.matrix[i][j] = sum;
      }
    return result;
  }      

  public double contraction()
  {
    int i,j;
    double result = 0d;
    for(i=0; i<rows; i++)
      for(j=0; j<columns; j++)
        result += matrix[i][j];
    return result;
  }  
  
  /**
   *  Return a matrix as a String.
   */
  public String toString()
  {
    if ((rows<=0) || (columns<=0))
      return "[]";
    int i,j;
    DecimalFormat format = new DecimalFormat("00.0000");
    format.setPositivePrefix("+");

    StringBuffer str = new StringBuffer();
    for(i=0; i<(rows-1); i++)
    {
      for(j=0; j<(columns-1); j++)
        if (Math.round(matrix[i][j]*10000)!=0)
          str.append(format.format(matrix[i][j])+" ");
        else
          str.append("-------- ");
      if (Math.round(matrix[i][columns-1]*10000)!=0)
        str.append(format.format(matrix[i][columns-1])+"\n");
      else
          str.append("--------\n");
    }
    for(j=0; j<(columns-1); j++)
      if (Math.round(matrix[rows-1][j]*10000)!=0)
        str.append(format.format(matrix[rows-1][j])+" ");
      else
          str.append("-------- ");
    if (Math.round(matrix[rows-1][columns-1]*10000)!=0)
      str.append(format.format(matrix[rows-1][columns-1]));
    else
      str.append("-------- ");
    return str.toString();
  }


  /** 
   * Diagonalize this matrix with the Jacobi algorithm.
   *
   * @param nrot Count of max. rotations
   * @return Matrix m, with m^t * this * m = diagonal
   *
   * @cdk.keyword Jacobi algorithm
   * @cdk.keyword diagonalization
   */
  public Matrix diagonalize(int nrot)
  {
    Matrix m = duplicate();
    if (m.rows!=m.columns)
        
    {
      System.err.println("Matrix.diagonal: Sizes mismatched");
      return null;
    }
    int n = m.rows;

    int j,iq,ip,i;
    
    double tresh,theta,tau,t,sm,s,h,g,c;
    double[] b,z;

    Matrix v = new Matrix(columns,columns);
    Vector d = new Vector(columns);
    
    b = new double[n+1];
    z = new double[n+1];
    for (ip=0;ip<n;ip++)
    {
      for (iq=0;iq<n;iq++)
        v.matrix[ip][iq]=0.0;
      v.matrix[ip][ip]=1.0;
    } 
    
    for (ip=0;ip<n;ip++)
    {
      b[ip]=d.vector[ip]=m.matrix[ip][ip];
      z[ip]=0.0;
    }

    nrot=0;
    for (i=1;i<=50;i++)
    {
      sm=0.0;
      for (ip=0;ip<n-1;ip++)
      {
        for (iq=ip+1;iq<n;iq++)
          sm += Math.abs(m.matrix[ip][iq]);
      }   

      // Ready ??
      if (sm == 0.0)
        return v;

      if (i < 4)
        tresh=0.2*sm/(n*n);
      else
        tresh=0.0;

      for (ip=0;ip<n-1;ip++)
      {
        for (iq=ip+1;iq<n;iq++)
        {
          g=100.0*Math.abs(m.matrix[ip][iq]);
          if ((i > 4) && (Math.abs(d.vector[ip])+g == Math.abs(d.vector[ip]))
                      && (Math.abs(d.vector[iq])+g == Math.abs(d.vector[iq])))
            m.matrix[ip][iq]=0.0;
          else if (Math.abs(m.matrix[ip][iq]) > tresh)
          {
            h = d.vector[iq]-d.vector[ip];
            if (Math.abs(h)+g == Math.abs(h))
              t = (m.matrix[ip][iq])/h;
            else 
            {
              theta = 0.5*h/(m.matrix[ip][iq]);
              t = 1.0/(Math.abs(theta)+Math.sqrt(1.0+theta*theta));
              if (theta < 0.0) t = -t;
            } 
            c = 1.0/Math.sqrt(1+t*t);
            s = t*c;
            tau = s/(1.0+c);
            h = t*m.matrix[ip][iq];
            z[ip] -= h;
            z[iq] += h;
            d.vector[ip] -= h;
            d.vector[iq] += h;
            m.matrix[ip][iq]=0.0;

            // Case of rotaions 1<=j<p
            for (j=0;j<ip;j++)
            {
              g=m.matrix[j][ip];
              h=m.matrix[j][iq];
              m.matrix[j][ip]=g-s*(h+g*tau);
              m.matrix[j][iq]=h+s*(g-h*tau);
            } 
            // Case of rotaions p<j<q
            for (j=ip+1;j<iq;j++) 
            {
              g=m.matrix[ip][j];
              h=m.matrix[j][iq];
              m.matrix[ip][j]=g-s*(h+g*tau);
              m.matrix[j][iq]=h+s*(g-h*tau);
            } 
            // Case of rotaions q<j<=n
            for (j=iq+1;j<n;j++) 
            {
              g=m.matrix[ip][j];
              h=m.matrix[iq][j];
              m.matrix[ip][j]=g-s*(h+g*tau);
              m.matrix[iq][j]=h+s*(g-h*tau);
            } 


            for (j=0;j<n;j++)
            {
              g=v.matrix[j][ip];
              h=v.matrix[j][iq];
              v.matrix[j][ip]=g-s*(h+g*tau);
              v.matrix[j][iq]=h+s*(g-h*tau);
            } 
            ++nrot;
          } 
        } 
      } 

      for (ip=0;ip<n;ip++)
      {
        b[ip] += z[ip];
        d.vector[ip]=b[ip];
        z[ip]=0.0;
      } 
    } 
    System.out.println("Too many iterations in routine JACOBI");
    return v;
  }

  /**
   * Solves a linear equation system with Gauss elimination.
   *
   * @cdk.keyword Gauss elimination
   */
  public static Vector elimination(Matrix matrix, Vector vector)
  {
    int i,j,k,ipvt;
    int n = vector.size;
    int[] pivot = new int[n];
    double c, temp;
    //double[] x = new double[n];
    Vector result = new Vector(n);
    Matrix a = matrix.duplicate();
    Vector b = vector.duplicate();
    
    for(j=0; j<(n-1); j++)
    {
      c = Math.abs(a.matrix[j][j]);
      pivot[j] = j;
      ipvt = j;
      for(i=j+1; i<n; i++)
        if (Math.abs(a.matrix[i][j])>c)
        {
          c = Math.abs(a.matrix[i][j]);
          ipvt = i;
        } 
        
      // Exchanges rows when necessary
      if (pivot[j]!=ipvt)
      {
        pivot[j] = ipvt;
        pivot[ipvt] = j;
        for(k=0; k<n; k++)
        {
          temp = a.matrix[j][k];
          a.matrix[j][k] = a.matrix[pivot[j]][k];
          a.matrix[pivot[j]][k] = temp;
        } 
        
        temp = b.vector[j];
        b.vector[j] = b.vector[pivot[j]];
        b.vector[pivot[j]] = temp;
      } 
      
      // Store multipliers
      for(i=j+1; i<n; i++)
        a.matrix[i][j] = a.matrix[i][j] / a.matrix[j][j];
        
      // Give elements below the diagonal a zero value
      for(i=j+1; i<n; i++)
      {
        for(k=j+1; k<n; k++)
          a.matrix[i][k] = a.matrix[i][k] - a.matrix[i][j]*a.matrix[j][k];
        b.vector[i] = b.vector[i] - a.matrix[i][j]*b.vector[j];
        
        a.matrix[i][j] = 0; // Not necessary
      } 
    } 
    
    // Rueckwaertseinsetzen (which is?)
    result.vector[n-1] = b.vector[n-1]/a.matrix[n-1][n-1];
    for(j=n-2; j>=0; j--)
    {
      result.vector[j] = b.vector[j];
      for(k=n-1; k>j; k--)
        result.vector[j] = result.vector[j] - result.vector[k]*a.matrix[j][k];
      result.vector[j] = result.vector[j]/a.matrix[j][j];
    } 
    
    return result;
  }

  /**
   * Orthonormalize the vectors of this matrix by Gram-Schmidt.
   *
   * @cdk.keyword orthonormalization
   * @cdk.keyword Gram-Schmidt algorithm
   */
  public Matrix orthonormalize(Matrix S)
  {
    int p,q,k,i,j;
    double innersum;
    double length;
    //Matrix scr = S.mul(this);
    Matrix result = duplicate();
    for(p=0; p<columns; p++) // Loops over all vectors
    {
      for(i=0; i<rows; i++)
        result.matrix[i][p] = matrix[i][p];

      for(k=0; k<p; k++)  // Substracts the previous vector 
      {
        // First the calculation of the product <phi_p|phi_k>=length
        length = 0;
        for(i=0; i<rows; i++) // Loops over all vectors
        { 
          innersum = 0;
          for(j=0; j<rows; j++)
          {
            innersum += result.matrix[j][p]*S.matrix[i][j];
          }
          length += result.matrix[i][k]*innersum;
        } 

        // Then the substraction of  phi_k*length
        for(q=0; q<rows; q++)
          result.matrix[q][p] -= result.matrix[q][k]*length;
      }

      // Calculates the integral for normalization
      length = 0;
      for(i=0; i<rows; i++)
        for(j=0; j<rows; j++)
          length += result.matrix[i][p]*result.matrix[j][p]*S.matrix[i][j];

      length = Math.sqrt(length);

      // Normalizes the vector
      if (length!=0d)
        for(q=0; q<rows; q++)
          result.matrix[q][p] /= length;
      else
        System.out.println("Warning(orthonormalize):"+(p+1)+". Vector has length null");
    }
    return result;
  }

  /**
   * Normalizes the vectors of this matrix.
   */
  public Matrix normalize(Matrix S)
  {
    int p,q,i,j;
    double length;
    Matrix result = duplicate();
    for(p=0; p<columns; p++) // Loops over all vectors
    {
      // Calculates the normalization factor
      length = 0;
      for(i=0; i<rows; i++)
        for(j=0; j<rows; j++)
          length += result.matrix[i][p]*result.matrix[j][p]*S.matrix[i][j];
          
      length = Math.sqrt(length);
      
      // Normalizes the vector
      if (length!=0d)
        for(q=0; q<rows; q++)
          result.matrix[q][p] /= length;
      else
        System.out.println("Warning(orthonormalize):"+(p+1)+". Vector has length null");
    }
    return result;
  }

}