#!/bin/sh
# the next line is executed by /bin/sh, but not tcl \
echo Tcl script name: $0, with $# arguments; echo ; exec tclsh "$0" "$@"

#*******************************************
# Call necessary pacakges
#*******************************************
package require math::linearalgebra

#*******************************************
# Definition of procedures
#*******************************************
proc det matrix {
    if {[llength $matrix] != [llength [lindex $matrix 0]]} {
       error "non-square matrix"
    }
    switch [llength $matrix] {
       2 {
           foreach {a b c d} [join $matrix] break
           expr {$a*$d - $b*$c}
       } default {
           set i 0
           set mat2 [lrange $matrix 1 end]
           set res 0
           foreach element [lindex $matrix 0] {
              if $element {
                 set sign [expr {$i%2? -1: 1}]
                 set res [expr {$res + $sign*$element* [det [cancelcol $i $mat2]]}]
              }
              incr i
           }
           set res
       }
    }
 }
 proc cancelcol {n matrix} {
    set res {}
    foreach row $matrix {
        lappend res [lreplace $row $n $n]
    }
    set res
 }



#*******************************************
# Main program starts here
#*******************************************

# Durer's magic square
puts "Durer's magic square"
set A {{16 3 2 13} {5 10 11 8} {9 6 7 12} {4 15 14 1}}
puts "A = \n[::math::linearalgebra::show $A %2d]"

# Transpose the matrix A
puts "Transpose of the matrix A, At ="
set At [::math::linearalgebra::transpose $A]
puts "[::math::linearalgebra::show $At %2d]"

# Matrix Summation
puts "A + At ="
puts "[::math::linearalgebra::show \
      [::math::linearalgebra::add $A $At] %2d]"

# Matrix multiplication
puts "At * A ="
set AtA [::math::linearalgebra::matmul $At $A]
puts "[::math::linearalgebra::show $AtA ]"

# Is a given matrix symmetric?
puts "isSymmetric( At * A ) = [::math::linearalgebra::symmetric $AtA]"

# Matrix determinant
set detA [det $A]
puts "det(A) = $detA.\n"

# Create a 3X3 matrix O whose elements are all zeros. 
puts "Create a 3X3 matrix O whose elements are all zeros."
set O [::math::linearalgebra::mkMatrix 3 3 0]
puts "O = \n[::math::linearalgebra::show $O]"

# Create a 3X3 identity matrix I
puts "Create a 3X3 identity matrix I."
set I [::math::linearalgebra::mkIdentity 3]
puts "I = \n[::math::linearalgebra::show $I]"

# Create a 3X3 random matrix R
puts "Create a 3X3 random matrix R."
set R [::math::linearalgebra::mkRandom 3]
puts "R = \n[::math::linearalgebra::show $R]"

# Get the 2nd column from the matrix R.
puts "Get the 2nd column from the matrix R."
set c2 [::math::linearalgebra::getcol $R 1]
puts "c2 = \n[::math::linearalgebra::show $c2]"