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
