namespace import ::math::linearalgebra::*

 proc trace {arg} {if {0} {puts "$arg"}}
 proc error {arg} {if {1} {puts "$arg"}}

 proc simplex_phase2 {l} {
     ## prepare matrix (fix float numbers)
     lassign $l cVect cMatrix bVect
     lassign [shape $cMatrix] rows cols
     set nvars [llength $cVect]
     set nconstr [llength $bVect]
     set l [list [scale_vect 1.0 $cVect] [scale_mat 1.0 $cMatrix] [scale_vect 1.0 $bVect]]
     trace "cVect:\n[show $cVect]\ncMatrix:\n[show $cMatrix]\nbVect:\n[show $bVect]"
     if {![expr {$nvars==$cols&&$nconstr==$rows}]} {
         error "problem is malformed (unmatching rows/cols count)."
         exit
     }


     set fix_b_sign 0
     for {set row 0} {$row < $rows} {incr row} {
         set bVect_i [lindex $bVect $row]
         if {$bVect_i < 0} {
             set cMatrix_i [getrow $cMatrix $row]
             lset cMatrix $row [scale -1.0 $cMatrix_i]
             lset bVect $row [expr -1.0*$bVect_i]
             set fix_b_sign 1
         }
     }

     if $fix_b_sign {trace "fixed cMatrix/bVect:\n[show $cMatrix]\n[show $bVect]"}
     
     ## find identity matrix (if any), that is:
     ## choose indices for in-base vars, and out-base vars
     lassign [shape $cMatrix] nc n
     set I [mkIdentity $nc]
     set iB {}
     set iN {}
     set cMt [transpose $cMatrix]
     foreach Ir $I {
         set idx [lsearch -exact $cMt $Ir]
         if {$idx != -1} {lappend iB $idx}
     }
     for {set iout 0} {$iout < $n} {incr iout} {
         if {[lsearch -exact $iB $iout] == -1} {
             lappend iN $iout
         }
     }
     unset cMt I Ir

     trace "iB=$iB, iN=$iN"

     if {[llength $iB] < $rows} {
         trace "### not ready for startup; entering phase 1"
         set new_cVect {}
         set n_artificial_vars 0
         set new_cMatrix {}
         set cMt [transpose $cMatrix]
         set cnt 0
         foreach cMatrix_i $cMatrix {
             set e_k {}
             for {set j 0} {$j < [llength $cMatrix]} {incr j} {
                 lappend e_k [expr $j==$cnt?1.0:0.0]
             }
             incr cnt
             if {[lsearch -exact $cMt $e_k] == -1} {
                 incr n_artificial_vars
                 lappend cMt $e_k
             }
         }
         set new_cMatrix [transpose $cMt]
         for {set j 0} {$j < $cols} {incr j} {
             lappend new_cVect 0.0
         }
         for {set j 0} {$j < $n_artificial_vars} {incr j} {
             lappend new_cVect 1.0
         }
         set result [simplex_phase2 [list $new_cVect $new_cMatrix $bVect]]
         array set _ [set result]
         trace "### end phase1: $result"
         for {set j $cols} {$j < [expr $cols+$n_artificial_vars]} {incr j} {
             if [expr [lindex $_(result:) $j]!=0.0?1:0] {
                 #TODO: continue swapping artificial vars with original vars
                 trace "### unfeasible problem - artificial vars remain in base"
                 return -1
             }
         }
         return $result
     }

     trace "problem ready for startup"

     ## algorithm startup
     set Bt {}
     set cb {}
     foreach ib $iB {
         lappend Bt [getcol $cMatrix $ib]
         lappend cb [lindex $cVect $ib]
     }
     set B [transpose $Bt]
     # inversion not required for first step, cause we have an identity matrix
     #set Binv [invert $B]
     set Binv $B
     set Nt {}
     set cn {}
     foreach in $iN {
         lappend Nt [getcol $cMatrix $in]
         lappend cn [lindex $cVect $in]
     }
     set N [transpose $Nt]
     set BinvN [matmul $Binv $N]
     set Binvb [matmul $Binv $bVect]

     set iteration 0
     set iB_history {}
     ### LOOP ENTER POINT
     while {1} {
         trace "*********** Iteration #[incr iteration] ******************"

         ## compute gamma (reduced costs) vector
         set gamma [axpy -1 [matmul $cb $BinvN] $cn]
         trace "gamma:\n[show $gamma]"

         set h 0
         foreach gamma_i $gamma {if {$gamma_i < 0} {break} {incr h}}
         ## check for optimal solution
         if {$h >= [llength $gamma]} {
             ## we found an optimal solution
             ## check uniqueness:
             set unique 1
             foreach gamma_i $gamma {if {$gamma_i == 0} {set unique 0}}

             ## now let's build our x vector (x_star)
             set x_star {}
             for {set ix 0} {$ix < $cols} {incr ix} {
                 set ib_idx [lsearch -exact $iB $ix]
                 if {$ib_idx != -1} {
                     lappend x_star [lindex $Binvb $ib_idx]
                 } else {
                     lappend x_star 0.0
                 }
             }
             #warn "found optimal solution: $x_star"
             #warn "uniqueness: $unique"
             #warn "value at optimum: [dotproduct $x_star $cVect]"
             return [list result: $x_star iB: $iB gamma: $gamma]
         } else {
             ## or unbound problem:
             for {set hi 0} {$hi < [llength $gamma]} {incr hi} {
                 ## foreach negative component of gamma, check if the B^(-1)N matching
                 ## column is all < 0
                 if {[lindex $gamma $hi] < 0} {
                     set BinvN_h [getcol $BinvN $hi]
                     set BinvN_neg_cnt 0
                     foreach BinvN_h_i $BinvN_h {if {$BinvN_h_i < 0} {incr BinvN_neg_cnt}}

                     if {$BinvN_neg_cnt == [llength $BinvN_h]} {
                         #warn "found unbound problem"
                         return -1
                     }
                 }
             }
         }

         ## - ingoing var and outgoing var -
         ## find h and k (using Bland's anti-loop rule)
         set hmin [lindex $gamma 0]
         set h 0
         for {set j 1} {$j < [llength $gamma]} {incr j} {
             set newhmin [lindex $gamma $j]
             if {$newhmin < [lindex $hmin 0]} {
                 set h $j
                 set hmin $newhmin
             } elseif {$newhmin < [lindex $hmin 0]} {
                 lappend h $j
                 lappend hmin $newhmin
             }
         }
         set h [lindex $h 0]
         set hmin [lindex $hmin 0]
         
         set pi_h [getcol $BinvN $h]
         
         set kmin [expr [lindex $Binvb 0]/[lindex $pi_h 0]]
         set k 0
         for {set j 1} {$j < [llength $pi_h]} {incr j} {
             set newkmin [expr [lindex $Binvb $j]/[lindex $pi_h $j]]
             if {($newkmin < $kmin && $newkmin >= 0) || $kmin < 0} {
                 set k $j
                 set kmin $newkmin
             }
         }

         trace "h=$h, k=$k, rho=$kmin"

         ## compute new iB,iN (indices for in-base/out-base vars)
         ## swapping h-th in-base var with k-th out-base var
         set iB_old $iB
         set iN_old $iN
         lset iB $k [lindex $iN_old $h]
         lset iN $h [lindex $iB_old $k]

         trace "iB=$iB, iN=$iN"

         ## pivot operation - prepare M matrix

         set M {}
         lappend M [getcol $BinvN $h]
         lassign [shape $BinvN] BinvN_rows BinvN_cols
         for {set ij 0} {$ij < $BinvN_cols} {incr ij} {
             if {$ij == $h} {
                 ## e_k here
                 set e_k [getcol $BinvN $h]
                 set e_k_i 0
                 foreach e_k___ $e_k {
                     lset e_k $e_k_i [expr ($e_k_i==$k)?1.0:0.0]
                     incr e_k_i
                 }
                 lappend M $e_k
             } else {
                 set BinvN_col [getcol $BinvN $ij]
                 lappend M $BinvN_col
             }
         }
         lappend M $Binvb
         set M [transpose $M]

         trace "before pivot:\nM:\n[show $M]"

         ## pivot operation:
         set row_k [getrow $M $k]
         lset M $k [scale [expr 1/[lindex [getcol $M 0] $k]] $row_k]
         set row_k_neg [scale -1.0 [getrow $M $k]]
         for {set ij 0} {$ij < [llength $M]} {incr ij} {
             if {$ij == $k} {continue}
             set cur_row [getrow $M $ij]
             lset M $ij [add_vect $cur_row [scale [lindex $cur_row 0] $row_k_neg]]
         }

         trace "after pivot:\nM:\n[show $M]"

         ## new BinvN, Binvb from M-pivoted
         set BinvN {}
         for {set ij 1} {$ij <= $BinvN_cols} {incr ij} {
             lappend BinvN [getcol $M $ij]
         }
         set BinvN [transpose $BinvN]
         set Binvb [getcol $M end]

         set cb {}
         set cn {}
         set Bt {}
         set Nt {}
         foreach ib $iB {
             lappend Bt [getcol $cMatrix $ib]
             lappend cb [lindex $cVect $ib]
         }
         foreach in $iN {
             lappend Nt [getcol $cMatrix $in]
             lappend cn [lindex $cVect $in]
         }
         set B [transpose $Bt]
         set N [transpose $Nt]

         trace "new BinvN:\n[show $BinvN]"
         trace "new Binvb:\n[show $Binvb]"
     } 
 }

 set f {-4 -1 -1}
 set c {
  {2 1 2}
  {3 3 1}
 }
 set b {4 3}
 set l [list $f $c $b]

 set result [simplex_phase2 $l]

 if {"$result" == "-1"} {
     puts "unfeasible or unbound problem"
 } else {
     puts "$result"
 }