#  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.DTRSM
import java.lang.Runtime

module ParallelMatrix

    PARALLEL_THRESHOLD = 10
    @@threads = 1
    
    @@parallel_multiply = false
      
        
    # Allows for the creation of class methods
    module ClassMethods; end
    def self.included(klass)
        klass.extend(ClassMethods)
    end
    
    def pmultiply_1d_byrow(m,threads)
        thr_hold = []
        rows = Array.new(self.row_size,nil)
        
        self.row_size.times do |i|

            completed = finished_thread(thr_hold,threads)
            if not completed.nil?
                rows[completed['rid']] = completed.value
                thr_hold.delete(completed)
            end
                
            thr_hold << Thread.new(i) {
                row = []
            
                m.column_size.times do |j|
                    row << self.row(i).dot(m.column(j))
                end
                
                row
            }
            thr_hold[-1]['rid'] = i
        end
        
        thr_hold.each do |t|
            t.join
            rows[t['rid']] = t.value
        end
        
        Matrix.rows(rows)   
    end

    def pmultiply_1d_bycolumn(m,threads)
        thr_hold = []
        columns = Array.new(m.column_size,nil)
        
        m.column_size.times do |j|

            completed = finished_thread(thr_hold,threads)
            if not completed.nil?
                columns[completed['rid']] = completed.value
                thr_hold.delete(completed)
            end
                
            thr_hold << Thread.new(j) {
                column = []
            
                self.row_size.times do |i|
                    column << self.row(i).dot(m.column(j))
                end
                
                column
            }
            thr_hold[-1]['rid'] = j
        end
        
        thr_hold.each do |t|
            t.join
            columns[t['rid']] = t.value
        end
        
        Matrix.columns(columns)        
    end
    
    def pmultiply_1d(m,threads)
        if(self.row_size > threads and self.row_size <= m.col_size)
            return pmultiply_1d_byrow(m,threads)
        elif(m.col_size > threads)
            return pmultiply_1d_bycolumn(m,threads)
        else
            return pmultiply_1d_byrow(m,threads)
        end        
    end
    
    def *(m)
    	  if(Matrix.is_parallel_multiply? and m.is_a? Matrix and @@threads > 1)
           	if(self.row_size > PARALLEL_THRESHOLD or m.column_size > PARALLEL_THRESHOLD)
               	pmultiply(m,@@threads)
           	else
               	multiply_single(m)
           	end
        else
           	multiply_single(m)
        end
    end

    def pinverse(threads)
        n = self.row_size
        raise ArgumentError if self.column_size != n
        
        # First generate the powers of the matrix
        # This may be heavy on the memory...
        powers = []
        powers << Matrix.identity(self.row_size)
        n.times do
            powers << powers[-1]*self
        end
        
        # Next compute all traces
        traces_mat = Matrix.diagonal(*(1...n)).to_java
        traces_vec = Matrix.zeros(n,1).to_java
        (1...n).each do |i|
            traces_vec[i][0] =  powers[i].trace
        end
        (1..n).each do |i|
            (n-i).times do |j|
                traces_mat[i+j][n-(i+j)] = traces_vec[i][0]
            end
        end
        
        # Retrieve the coefficients of the characteristic polynomial
        # via a back-subsitution
        org.netlib.blas.DTRSM('L','L','N','N',n,1,1.0,traces_mat,traces_vec)
        
        # Sum the series to retrieve the inverse
        ainv = powers[n-1]
        (1..n).each do |i|
            ainv = ainv + traces_vec[i][0]*powers[n-i-1]
        end
        
    end

    def finished_thread(thread_list, thread_count)
        return nil if(thread_list.size < thread_count)
        
        alive = true
        while(alive) do
            thread_list.each do |t|
                return t unless t.alive?
            end
        end
    end

    module ClassMethods

        def set_thread_count(*args)
            if(args.nil?)
              @@threads = java.lang.Runtime.getRuntime().availableProcessors()
            elsif(args.length != 1)
              @@threads = java.lang.Runtime.getRuntime().availableProcessors()
            elsif(args.length == 1)
              @@threads = args[0] if args[0].is_a? Fixnum
            else
              @@threads = 1
            end
            @@threads = 1 unless @@threads > 0
        end
        
        def get_thread_count
            return @@threads          
        end
        
        def set_parallel_multiply(value)
            @@parallel_multiply = value
        end
        
        def is_parallel_multiply?
            return @@parallel_multiply
        end
    end
    
end

class Matrix
    alias_method :multiply_single, :*
    include ParallelMatrix
end  
