Advanced Search

Venkat Subramanian

386 Reputation

13 Badges

15 years, 240 days
Contact 119368

MaplePrimes Activity


These are questions asked by

Calling UMFPACK sparsedirect linearsolve...

Dr. Venkat Subramanian 386 Maple
sparse linear-algebra
February 16 2022
0 2

@acer,
In the past, you had shown how to use take a given sparse A matrix and store the LU components and use the same factorization for multiple backsolves at https://www.mapleprimes.com/posts/41191-Solving-Sparse-Linear-Systems-In-Maple#comment200817

When PDEs are solved, we are calling the A matrix at every time step (with some discretizations in x) and the factorization is done and stored at every time step (this is time-consuming and memory-consuming). Is it possible to call UMFPACK using only the sparse storage and entries? The main routine seems to provide the option 

https://people.sc.fsu.edu/~jburkardt/f77_src/umfpack/umfpack.html

This would mean that the pattern is found only once for the Jacobian at t= 0, then only the non-zero sparse entries (vector/row, not a matrix) are updated at every time step.

To be clear, what I am asking for is create a random sparse matrix (say 4x4 or 10x10)
create a b Vector.

Solve AZ =b with method = SparseDirect. (This should internally create and store R, CC, X as at https://people.sc.fsu.edu/~jburkardt/f77_src/umfpack/umfpack.html)

Store R, CC, X and update only X for a different A matrix with the same sparsity pattern and solve for Z using stored (R, CC), and updated X.

Non symmetric sparse solver - non-implem...

Dr. Venkat Subramanian 386 Maple
sparse linearalgebra linearsolve
May 06 2020
2 1

Hi All,

https://www.maplesoft.com/support/helpJP/Maple/view.aspx?path=NAG%2ff11dec

This link shows that there are 3 options rgmres, cgs and bicgstab. When I use iterative solver, it uses the cgs option.

LinearAlgebra:-LinearSolve(
      A, b, method= SparseIterative)

LinearSolve:   "copying right hand side to enable external call"
LinearSolve:   "using method"   SparseIterative
LinearSolve:   "using method  "   SparseIterative
LinearSolve:   "calling external function"
LinearSolve:   "using CGS method"
LinearSolve:   "preconditioning with incomplete LU factorization"
LinearSolve:   "level of fill = "   0
LinearSolve:   "using complete pivoting strategy"
LinearSolve:   "dimension of workspaces for preconditioner = "   44
LinearSolve:   "using infinity norm in stopping criteria"
LinearSolve:   "setting maximum iterations to "   200
LinearSolve:   "setting tolerance to "   .10e-7
LinearSolve:   "NAG"   hw_f11zaf
LinearSolve:   "NAG"   hw_f11daf
LinearSolve:   "NAG"   hw_f11dcf
LinearSolve:   "number of iterations"   1
LinearSolve:   "residual computed last as"   3.33066907387547e-016
 

 

How can I force Maple to use BiCGSTAB?

thanks

Is it possible to edit or modify mla fil...

Dr. Venkat Subramanian 386 Maple
debugger printlevel trace
December 08 2016
1 3

I understand how Browse() works. Sometimes Maple crashes or fails for some programs. Is it possible to add print("here" ) into the the mla files in the library to identify why and where Maple fails?

 

Thanks

 

Recovering un-returned local variables f...

Dr. Venkat Subramanian 386 Maple
procedure local showstat
December 08 2016
1 5

I would like to return local variable y (line 4 in showstat) in the attached dummy procedure (s1) without manually adding any comment inside the procedure s1. This procedure is a simple one and easy to copy paste/or change. When we have a long procedure, it is difficult to do so. I will always know the name of the local variable I want (say, y) and/or line number in showstat

Thanks

PS, I want to get y:=array(1..,2[(1)=x,2=zz])
 

> 

restart;
> 

s1:=proc(n,x)
local y,xx,i,j,zz::array(1..n,1..n);
for i from 1 to n do for j from 1 to n do zz[i,j]:=x[i]*(1+x[j]^2);od:
od:
y:=array(1..2,[(1)=x, (2)=zz]):
for j from 1 to n do xx[i]:=zz[i,i]/(add(zz[i,j],j=1..n));od:
0;
end proc;

s1 := proc (n, x) local y, xx, i, j, zz::(array(1 .. n, 1 .. n)); for i to n do for j to n do zz[i, j] := x[i]*(1+x[j]^2) end do end do; y := array(1 .. 2, [1 = x, 2 = zz]); for j to n do xx[i] := zz[i, i]/add(zz[i, j], j = 1 .. n) end do; 0 end proc

 (1)
> 

showstat(s1);

 

s1 := proc(n, x)

local y, xx, i, j, zz::array(1 .. n,1 .. n);

   1   for i to n do

   2     for j to n do

   3       zz[i,j] := x[i]*(1+x[j]^2)

         end do

       end do;

   4   y := array(1 .. 2,[1 = x, 2 = zz]);

   5   for j to n do

   6     xx[i] := zz[i,i]/add(zz[i,j],j = 1 .. n)

       end do;

   7   0

end proc

 

 
> 

x0:=Vector(2,[1,1]);

x0 := Vector(2, {(1) = 1, (2) = 1})

 (2)
> 

s1(2,x0);

0

 (3)
> 

 


 

Download showstatexample.mws

Bug in dsolve numeric or fdiff?...

Dr. Venkat Subramanian 386
dsolve numeric differential-equations
November 03 2016
1 2


To summarize,

dsolve numeric parametric form gives wrong answers when fdiff is used to calculate the Hessian. Doing dsolve twice gives the correct answer. Hope I am not making syntax or programming errors.

Also, fdiff is not compatible with vector form.

> 

restart;
> 

Digits:=15;

Digits := 15

 (1)
> 

sys:=diff(y1(t),t)=-(u+u^2/2)*y1(t),diff(y2(t),t)=u*y1(t);

sys := diff(y1(t), t) = -(u+(1/2)*(u^2))*y1(t), diff(y2(t), t) = u*y1(t)

 (2)

 The system is solved with initial conditions 1,0 and value of u being equal to u1 for t<=0.5, and u2 for t>0.5. The objective is y2 at t =1. There are errors when Hessian is calculated using parametric dsolve and fdiff

> 

sol1 := dsolve({sys, y1(0) = alpha, y2(0) = beta}, type = numeric, 'parameters' = [alpha, beta, u],maxfun=0,range=0..0.5):
> 

obj1:=proc(u1,u2)
local z1,s1,s2,z2;
global sol1;
sol1('parameters'=[1,0,u1]);
z1 := sol1(.5);
s1:=subs(z1,y1(t));
s2:=subs(z1,y2(t));
sol1('parameters'=[s1,s2,u2]):
z2:=sol1(0.5);
-subs(z2,y2(t));
end proc;

obj1 := proc (u1, u2) local z1, s1, s2, z2; global sol1; sol1('parameters' = [1, 0, u1]); z1 := sol1(.5); s1 := subs(z1, y1(t)); s2 := subs(z1, y2(t)); sol1('parameters' = [s1, s2, u2]); z2 := sol1(.5); -subs(z2, y2(t)) end proc

 (3)
> 

u0:=[0.8,1.8];

u0 := [.8, 1.8]

 (4)
> 

obj1(op(u0));

-.552540796143903

 (5)
> 

Hess1:=Matrix(2,2):for i from 1 to 2 do for j from 1 to 2 do Hess1[i,j]:=fdiff(obj1,[i,j],u0,workprec=1.0);od;od;Hess1;

Matrix(2, 2, {(1, 1) = 115.470000000000, (1, 2) = -1.11000000000000, (2, 1) = -1.11000000000000, (2, 2) = 2.22000000000000})

 (6)

 Hessian is wrong. It should be noted that when dsolve is called twice inside the objective function there is no error in Hessian calculation with fdiff. The parametric dsolve is efficient, how to get correct answers with that? Another question is how to use fdiff with vector form instead u1,u2 etc to make it generic for any number of variables?

> 

obj2:=proc(u1,u2)
local z1,s1,s2,z2,sol2;
sol2:=dsolve({op(subs(u=u1,[sys])), y1(0) = 1, y2(0) = 0}, type = numeric,maxfun=0,range=0..0.5);
z1 := sol2(.5);
s1:=subs(z1,y1(t));
s2:=subs(z1,y2(t));
sol2:=dsolve({op(subs(u=u2,[sys])), y1(0) = s1, y2(0) = s2}, type = numeric,maxfun=0,range=0..0.5):
z2:=sol2(0.5);
-subs(z2,y2(t));
end proc;

obj2 := proc (u1, u2) local z1, s1, s2, z2, sol2; sol2 := dsolve({op(subs(u = u1, [sys])), y1(0) = 1, y2(0) = 0}, type = numeric, maxfun = 0, range = 0 .. .5); z1 := sol2(.5); s1 := subs(z1, y1(t)); s2 := subs(z1, y2(t)); sol2 := dsolve({op(subs(u = u2, [sys])), y1(0) = s1, y2(0) = s2}, type = numeric, maxfun = 0, range = 0 .. .5); z2 := sol2(.5); -subs(z2, y2(t)) end proc

 (7)
> 

obj2(op(u0));

-.552540796143903

 (8)
> 

Hess2:=Matrix(2,2):for i from 1 to 2 do for j from 1 to 2 do Hess2[i,j]:=fdiff(obj2,[i,j],u0,workprec=1.0);od;od;Hess2;

Matrix(2, 2, {(1, 1) = .234295697631809, (1, 2) = 0.101868715004724e-1, (2, 1) = 0.101868715004724e-1, (2, 2) = 0.853480932660994e-1})

 (9)

 Doing dsolve twice gives correct answer compared to one parametric dsolve even though objective returns the same number (even gradient returns the expected answers).


 

Download Bugindsolveorfdiff.mws

1 2 3 Page 2 of 3