Actually I just noticed there is a fast solver in Lapack for solving Ax=b. This is all I need - thanks!
load 'math/lapack/gesv' m1=: ?.4 3$30 gesv_jlapack_ m;m1 NB. solves m * x = m1 Brian Babiak MD https://drbabiak.com/ <http://drbabiak.com> paypal: *http://www.paypal.me/BrianBabiakMD <https://www.paypal.me/BrianBabiakMD>* cash, check, debit/credit cards, paypal, apple pay, android pay, discover and amex accepted for all payments square: *https://squareup.com/store/brian-babiak-md <https://squareup.com/store/brian-babiak-md>* On Fri, Jul 28, 2017 at 10:03 AM, Brian Babiak <bdbab...@gmail.com> wrote: > Thanks for the help everyone. LAPACK is very useful. I also appreciate the > Rosetta code mesh. It's interesting how many choices were made the same. I > appreciate the help on this! > > Brian Babiak MD > https://drbabiak.com/ <http://drbabiak.com> > > paypal: *http://www.paypal.me/BrianBabiakMD > <https://www.paypal.me/BrianBabiakMD>* > > cash, check, debit/credit cards, paypal, apple pay, android pay, discover > and amex accepted for all payments > > square: *https://squareup.com/store/brian-babiak-md > <https://squareup.com/store/brian-babiak-md>* > > > On Fri, Jul 28, 2017 at 2:51 AM, Lippu Esa <esa.li...@varma.fi> wrote: > >> I think too. Here is the smallest eigenvalue of a: >> >> require 'math/lapack/geev' >> ev=.geev_jlapack_ a >> {.sort|>1{ev >> 6.83897e_14 >> >> Esa >> -----Original Message----- >> From: Programming [mailto:programming-boun...@forums.jsoftware.com] On >> Behalf Of 'Jon Hough' via Programming >> Sent: Friday, July 28, 2017 9:17 AM >> To: programm...@jsoftware.com >> Subject: Re: [Jprogramming] Matrix operations questions >> >> I think a is singular. You should use LAPACK to decompose it perhaps. >> see: http://code.jsoftware.com/wiki/Studio/LAPACK >> >> -------------------------------------------- >> On Fri, 7/28/17, Brian Babiak <bdbab...@gmail.com> wrote: >> >> Subject: Re: [Jprogramming] Matrix operations questions >> To: programm...@jsoftware.com >> Date: Friday, July 28, 2017, 2:33 PM >> >> Very helpful! At least I can use >> +/ .* for an iterative solution to Lp=I, solving for p. But >> why can’t I always get an inverse? (Even with non >> -infinite resistances, I get domain errors.) >> >> a=.L 2 2;3 1;1 4 >> a >> 4 _2 0 0 0 >> _2 0 0 0 0 0 0 0 0 >> 0 0 0 0 0 0 0 0 0 0 0 >> _2 8 _2 0 0 _2 _2 0 0 >> 0 0 0 0 0 0 0 0 >> 0 0 0 0 0 0 0 0 >> 0 _2 8 _2 >> 0 0 _2 _2 0 0 0 0 0 >> 0 0 0 0 0 0 0 0 0 0 0 >> 0 >> 0 0 _2 8 _2 0 0 _2 >> _2 0 0 0 0 0 0 0 0 >> 0 0 0 0 0 0 0 0 >> 0 0 0 >> _2 1004 0 0 0 _2 _1000 0 0 >> 0 0 0 0 0 0 0 0 0 0 0 >> 0 0 >> _2 _2 0 0 0 8 _2 0 >> 0 0 _2 0 0 0 0 0 >> 0 0 0 0 0 0 0 0 0 >> 0 _2 >> _2 0 0 _2 12 _2 0 0 _2 _2 >> 0 0 0 0 0 0 0 0 0 0 0 >> 0 0 >> 0 0 _2 _2 0 0 _2 1010 >> _2 0 0 _2 _1000 0 0 0 0 >> 0 0 0 0 0 0 0 0 >> 0 0 0 >> _2 _2 0 0 _2 2008 _1000 0 0 _1000 >> _2 0 0 0 0 0 0 0 0 0 0 >> 0 >> 0 0 0 0 _1000 0 0 0 >> _1000 4000 0 0 0 _1000 _1000 0 0 >> 0 0 0 0 0 0 0 0 >> 0 0 0 >> 0 0 _2 _2 0 0 0 8 _2 >> 0 0 0 _2 0 0 0 0 0 0 0 >> 0 0 >> 0 0 0 0 0 0 _2 >> _2 0 0 _2 1009 _1000 0 0 _2 >> _1 0 0 0 0 0 0 0 0 >> 0 >> 0 0 0 0 0 0 _1000 _1000 0 0 _1000 >> 5999 _1000 0 0 _999 _1000 0 0 0 0 0 0 >> 0 >> 0 0 0 0 0 0 0 0 >> _2 _1000 0 0 _1000 2008 _2 0 0 >> _2 _2 0 0 0 0 0 0 >> 0 0 0 >> 0 0 0 0 0 0 _1000 0 0 >> 0 _2 1006 0 0 0 _2 _2 0 0 0 0 >> 0 >> 0 0 0 0 0 0 0 0 >> 0 0 _2 _2 0 0 0 7 _1 >> 0 0 0 _2 0 0 0 0 >> 0 0 >> 0 0 0 0 0 0 0 0 0 _1 >> _999 0 0 _1 1004 _1 0 0 _1 _1 0 0 >> 0 >> 0 0 0 0 0 0 0 0 >> 0 0 0 0 _1000 _2 0 0 _1 >> 1009 _2 0 0 _2 _2 0 0 >> 0 0 0 >> 0 0 0 0 0 0 0 0 0 >> 0 _2 _2 0 0 _2 12 _2 0 0 _2 _2 >> 0 >> 0 0 0 0 0 0 0 0 >> 0 0 0 0 0 0 _2 0 >> 0 0 _2 8 0 0 0 _2 _2 >> 0 0 >> 0 0 0 0 0 0 0 0 0 0 >> 0 0 0 _2 _1 0 0 0 5 _2 0 >> 0 0 >> 0 0 0 0 0 0 0 >> 0 0 0 0 0 0 0 0 0 >> _1 _2 0 0 _2 7 _2 0 0 >> 0 >> 0 0 0 0 0 0 0 0 0 0 >> 0 0 0 0 0 0 _2 _2 0 0 _2 >> 8 _2 0 >> 0 0 0 0 0 0 0 >> 0 0 0 0 0 0 0 0 0 >> 0 0 _2 _2 0 0 _2 8 _2 >> 0 >> 0 0 0 0 0 0 0 0 0 0 >> 0 0 0 0 0 0 0 0 _2 0 >> 0 0 _2 4 >> -/ . * a >> _8.08478e26 >> %. a >> |domain error >> | %.a >> >> >> Brian Babiak >> MD >> https://drbabiak.com >> >> >> 167 Ridgecrest Rd. >> Ithaca >> NY 14850 – 9449 >> (607) 379 – 1335 >> fax (866) 813 – 4158 >> >> Payment links: >> >> http://PayPal.Me/BrianBabiakMD http://mkt.com/brian-babiak-md >> >> >> On 7/28/17, >> 1:03 AM, "Programming on behalf of 'Jon Hough' >> via Programming" <programming-boun...@forums.jsoftware.com >> on behalf of programm...@jsoftware.com> >> wrote: >> >> Matrix product >> is +/ .* (dyadic) >> determinant is -/ .* >> (monadic) >> >> see: >> http://code.jsoftware.com/wiki/Vocabulary/dot >> >> The matrix g is >> singular (since the rows are identical) and has no inverse. >> Hence %. gives an error. >> >> i.e. >> >> (-/ .*) g >> >> is zero. >> >> Hope that helps. >> >> >> Jon >> >> >> -------------------------------------------- >> On Fri, 7/28/17, Brian Babiak <bdbab...@gmail.com> >> wrote: >> >> >> Subject: [Jprogramming] Matrix operations questions >> To: programm...@jsoftware.com >> Date: Friday, July 28, 2017, 1:56 PM >> >> I’m working on >> a project to compute the >> distance >> resistance between two vertices in an electronic >> grid. I will share that later. But I am >> having domain error >> issues with >> matrix inverse and matrix product operations. >> See the following: >> >> >> >> >> f=.2 2$3 4 1 4 >> >> >> g=.2 2$i.2 >> >> >> h=.2 2$_3 4 2 _3 >> >> >> >> 2 >> 3$'f';'g';'h';f;g;h >> >> >> ┌───┬───┬─────┐ >> >> │f │g >> │h │ >> >> >> ├───┼───┼─────┤ >> >> │3 4│0 >> 1│_3 4│ >> >> >> │1 4│0 1│ 2 _3│ >> >> >> └───┴───┴─────� >> >> >> >> 2 >> >> 3$'$f';'$g';'$h';($f);($g);($h) >> >> >> ┌───┬───┬───┐ >> >> >> │$f │$g │$h │ >> >> >> ├───┼───┼───┤ >> >> >> │2 2│2 2│2 2│ >> >> >> └───┴───┴───� >> >> >> >> >> f . g >> >> |domain error >> >> | f .g >> >> >> f . h >> >> >> |domain error >> >> >> | f .h >> >> >> %. f >> >> >> 0.5 _0.5 >> >> >> _0.125 0.375 >> >> >> %. g >> >> >> |domain error >> >> >> | %.g >> >> >> f . f >> >> >> |domain error >> >> >> | f .f >> >> >> g . g >> >> >> |domain error >> >> >> | g .g >> >> >> h . h >> >> >> |domain error >> >> >> | h .h >> >> >> h * g >> >> >> 0 4 >> >> 0 _3 >> >> >> >> The answers to all these operations >> are >> simple. Why are there these >> domain errors? >> >> >> >> >> Here is >> the project this issue is >> coming up >> in, although I believe it is more basic than >> that. >> >> >> >> vw=: (,@(((_ >> >> 0)$~#@:])`]`(1$~,~@:}.@[)}))`(,@(((0 >> >> _)$~#@:])`]`(1$~,~@:}.@[)}))@.(<:@:({.@:[)) >> >> NB. (moveno >> boardsz) nw (x1 y1;x2 y2; >> etc) >> >> NB. Gives vertex >> weights for a move >> list of a game as >> vector >> >> >> >> >> >> rows=: >> <"1 >> >> >> cols=: <"1@|: >> >> diags=: (-.((0;_1)&{))@:(</.) >> >> a=: 3 3$i.9 >> >> conndiag=: 3 : >> 0 >> >> r=. r,.|.r >> [ r=. >0{z [ i=. 1 [ n=. >> #z [ z=. >> diags y >> >> >> while. i<n do. >> >> >> col=. >i{z >> >> i=. >:i >> >> adds=. >2<\col >> >> r=. r,adds >> >> r=. >> r,(|."1 adds) >> >> end. >> >> ) >> >> >> >> connsq=: 3 : 0 >> >> >> r=. 0 2$0 [ i=. 0 [ n=. #z [ z=. >> y >> >> while. >> i<n do. >> >> >> col=. >i{z >> >> >> i=. >:i >> >> >> adds=. >2<\col >> >> r=. r,adds >> >> r=. r,(|."1 adds) >> >> end. >> >> ) >> >> connrws=: >> connsq@:rows >> >> >> conncls=: connsq@:cols >> >> >> >> edges=: 3 : 0 NB. takes board size as >> argument >> >> a=. (2$bs)$i.*:bs [ bs=. y >> >> r=. (connrws >> a),(conncls a) >> >> >> r=. r,conndiag a >> >> >> ) >> >> >> >> BOARDSZ=: 5 >> >> EDGES=: edges >> BOARDSZ >> >> >> VPLYR=: 1 NB. whether the vertical >> >> player is first or second >> >> >> >> A=: 3 : 0 >> >> mat=. (2$netsz)$0 [ netsz=. *:BOARDSZ >> [ >> j=. 0 [ i=. 0 [ arg=. >> VPLYR,BOARDSZ >> >> >> if. y=0 do. pcweights=. netsz$1 else. >> >> pcweights=. arg vw y end. >> >> while. i<netsz do. >> >> >> while. j<netsz do. >> >> mat=. >> (0{(((1 2$i,j) e. >> >> EDGES)*((i{pcweights)+j{pcweights))) (<(i,j))}mat >> >> j=. >> >:j >> >> >> end. >> >> j=. >> 0 >> >> i=. >> >:i >> >> >> end. >> >> mat >> >> ) >> >> >> >> D=: 3 : 0 >> >> diag=. >> (2$netsz)$ 1,netsz#0 [ netsz=. >> >> *:BOARDSZ [ colsum=. +/y [ y=. A y >> >> >> r=. colsum*diag >> >> >> ) >> >> >> >> L=: D-A NB. the Laplacian >> >> >> >> NB. %. Laplacian >> or . (dot) for >> Laplacian times >> current gives domain errors >> >> >> >> Any help would be appreciated! >> >> >> >> Brian Babiak >> MD >> >> https://drbabiak.com >> >> >> >> >> >> >> >> ---------------------------------------------------------------------- >> For information about J forums see http://www.jsoftware.com/forum >> s.htm >> >> ---------------------------------------------------------------------- >> For information about J forums see http://www.jsoftware.com/forum >> s.htm >> >> >> ---------------------------------------------------------------------- >> For information about J forums see http://www.jsoftware.com/forums.htm >> ---------------------------------------------------------------------- >> For information about J forums see http://www.jsoftware.com/forums.htm >> ---------------------------------------------------------------------- >> For information about J forums see http://www.jsoftware.com/forums.htm >> > > ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
