#  JLMatrix - Matrix support for JRuby
#  Copyright (C) 2008 Approximatrix, LLC
#  http://www.approximatrix.com
#  
#  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 3 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, see <http://www.gnu.org/licenses/>.
#
#  NOTE: This software is licensed under the GPL and does not allow
#        linking to closed-source software.
#

import org.netlib.blas.DDOT
import org.netlib.blas.DAXPY
import org.netlib.blas.DNRM2
import org.netlib.blas.DGEMV
import org.netlib.blas.DGEMM
import org.netlib.lapack.DGETRF
import org.netlib.lapack.DGETRI
import org.netlib.lapack.DLANGE
import org.netlib.lapack.DGECON
import org.netlib.util.intW
import org.netlib.util.doubleW

require 'matrix'

include Math

class Matrix 
    alias_method :m, :row_size
    alias_method :n, :column_size

    # Converts the matrix to a 2D java array      
    def to_java
        ret = Java::double[m][n].new
        
        m.times do |i|
          n.times do |j|
            ret[i][j] = self[i,j]
          end
        end
        ret
    end
   
    alias_method :matrix_multiply, :*
      
    # Implements multiply using jLAPACK where appropriate      
    def *(x)
      if(x.is_a? Vector)
        raise ArgumentError, "Inner dimensions must match (#{m}x#{n} times #{x.length}x1)" unless x.length == n
        result = Vector.zeroes(m).to_java
        org.netlib.blas.DGEMV.DGEMV('N',m,n,1.0,to_java,x.to_java,1,0.0,result,1)
        return Vector.from_java result
      elsif(x.is_a? Matrix)
        raise ArgumentError, "Inner dimensions must match (#{m}x#{n} times #{x.m}x#{x.n})" unless x.m == n
        result = Java::double[m][x.n].new
        org.netlib.blas.DGEMM.DGEMM('N','N',m,x.n,n,1.0,to_java,x.to_java,1.0,result)
        return Matrix.from_java result 
      else
        matrix_multiply(x)
      end
        
    end

    # Raises all elements to a power
    def **(x)
        ret = []
        column_vectors.each do |ele|
            ret << (ele**x)
        end
        Matrix.columns(ret)
    end
    
    # Implements division by a scalar
    def /(value)
    	raise ArgumentError,"Only scalars may be divisors of Matrices" unless (value.is_a? Fixnum or value.is_a? Bignum or value.is_a? Float)
    	ret = []
    	column_vectors.each do |x|
    		ret << (x / value)
    	end
    	Matrix.columns(ret)
    end
  
    # Returns a Vector containing the diagonal elements of a symmetric matrix
    def diagonal(*args)
        raise NoMethodError, "Method only available for symmetric matrices" unless m==n
        holding = []
        m.times do |i|
            holding << self[i,i]
        end
        Vector.elements(holding)
    end
  
    # Computes the inverse using jLAPACK rather than the built-in Ruby routines
    def inverse
        raise NoMethodError, "Method only available for symmetric matrices" unless m==n
               
        temp = to_java
     
        ipiv = Java::int[m].new
        work = Java::double[m].new
        result = org.netlib.util.intW.new(0)
         
        # Together, these invert a matrix (temp)
        org.netlib.lapack.DGETRF.DGETRF(m,m,temp,ipiv,result);
        org.netlib.lapack.DGETRI.DGETRI(m,temp,ipiv,work,m,result);
             
        Matrix.from_java temp
    end
  
    # Computes the reciprocal of the matrix condition number
    def reciprocal_condition
        raise NoMethodError, "Method only available for symmetric matrices" unless m==n
  
        ipiv = Java::int[m].new
        work = Java::double[m].new
        one_norm_work = Java::double[m*4].new
        result = org.netlib.util.intW.new(0)
  
        # Determines an estimate of the reciprocal of the condition #
        temp_matrix = to_java
        one_norm = org.netlib.lapack.DLANGE.DLANGE("1",m,m,temp_matrix,work)
  
        temp_matrix = to_java
        org.netlib.lapack.DGETRF.DGETRF(m,m,temp_matrix,ipiv,result)
  
        temp_matrix = to_java
        rcond = org.netlib.util.doubleW.new(0.0);
        org.netlib.lapack.DGECON.DGECON("1",m,temp_matrix,one_norm,rcond,one_norm_work,ipiv,result)
  
        return 0.0 unless result.val == 0
        return rcond.val
  
     end
  
    alias_method :rcond, :reciprocal_condition
  
    # Constructs a Matrix from a Java 2D array
    def Matrix.from_java(x)
        rows = []
        begin
            x.length.times do |i|
                column = []
                x[i].length.times do |j|
                    column << x[i][j] 
                end
                rows << column
            end
            return Matrix.rows(rows)
        rescue
            raise ArgumentError,"Input must be a valid 2D Java array with consistent row lengths."
        end
        
    end
    
    # Creates a Matrix of zeros of arbitrary size
    def Matrix.zeroes(m,n)
    	rows = []
    	m.times do 
    		col = []
    		n.times do
    			col << 0.0
    		end
    		rows << col
    	end
    	Matrix.rows(rows)
    end     
        
end
