I'm using the thing for mesh resistance as someone noted above, The gesv seems to work very well, and a 121 x 121 resistor array was timed a 0.0008s (0.8 ms!) That's pretty amazing. %. doesn't seem to be anywhere close to that fast, and throws domain errors on many physical electrical network, whereas the gesv seems to be solid given the poking around I did a little with it today. Over the weekend, I'll have more time to test it out and solve. The code snippet I showed was for a Shannon Switching game, to just try out the Laplacian etc. I'll let you know how it works.
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 4:30 PM, Raul Miller <rauldmil...@gmail.com> wrote: > I can't test this because I do not know what you are using for m, and > because I get a value error on z2d when I try gesv_jlapack_ on > arbitrary values. > > Still, I am curious - if gesv_jlapack_ is better than (x %. y) what > are the conditions where it's better? > > Thanks, > > -- > Raul > > On Fri, Jul 28, 2017 at 1:00 PM, Brian Babiak <bdbab...@gmail.com> wrote: > > 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 > ---------------------------------------------------------------------- > For information about J forums see http://www.jsoftware.com/forums.htm > ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
