# Matrix-vector multiplication in weak formulation **URL:** \ **Category:** General Discussion\ **Created:** [May 31, 2020, 3:14pm UTC](https://community.freefem.org/t/matrix-vector-multiplication-in-weak-formulation/449 "2020-05-31T15:14:04Z")\ **Posts on this page:** 20\ **Page:** 1
**Author:** ![andryr](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/andryr/32/170_2.png) [@andryr](https://community.freefem.org/u/andryr)\ **Post date:** [May 31, 2020, 3:14pm UTC](https://community.freefem.org/t/matrix-vector-multiplication-in-weak-formulation/449/1 "2020-05-31T15:14:04Z")
Hello everyone I am having issues trying to solve a problem that has a matrix-vector product in its weak formulation. I defined a macro epsilon as follows : ```auto macro epsilon(u1, u2) [dx(u1), dy(u2), dy(u1)+dx(u2)] // ``` and there is a term in the weak formulation that looks like this : ```auto int2d(Th)( epsilon(uu2, vv2)' * C * epsilon(w, s) ) ``` where C is a 3x3 real matrix. The issue is that this code won’t compile and this is the error message I get ```auto 69 : (C * epsilon(uu2, vv2) [dx(uu2), dy( vv2), dy(uu2)+dx( vv2)] ) error operator * <14Matrice_CreuseIdE>, <10LinearCombI7MGauche4C_F0E> ``` However, when I use remove the matrix C from the weak formulation it compiles and run fine. ```auto int2d(Th)( epsilon(uu2, vv2)' * epsilon(w, s) ) ``` Any help would be appreciated Thanks ---
**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\ **Post date:** [May 31, 2020, 3:41pm UTC](https://community.freefem.org/t/matrix-vector-multiplication-in-weak-formulation/449/2 "2020-05-31T15:41:27Z")
Sadly, I don’t think it’s quite possible. However, you can trick the `varf` with a hand-defined macro. ```auto mesh Th; real[int, int] C(2,2); macro Cdense [[C(0,0), C(0,1)],[C(1,0), C(1,1)]]// varf vPb(u, v)= int2d(Th)( [2*v, v]'*Cdense*[3*u, u] ); ``` ---
**Author:** ![andryr](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/andryr/32/170_2.png) [@andryr](https://community.freefem.org/u/andryr)\ **Post date:** [May 31, 2020, 4:09pm UTC](https://community.freefem.org/t/matrix-vector-multiplication-in-weak-formulation/449/3 "2020-05-31T16:09:22Z")
It works ! Thank you very much. ---
**Author:** ![Lee](https://avatars.discourse-cdn.com/v4/letter/l/bc8723/32.png) [@Lee](https://community.freefem.org/u/Lee)\ **Post date:** [July 10, 2020, 3:44am UTC](https://community.freefem.org/t/matrix-vector-multiplication-in-weak-formulation/449/4 "2020-07-10T03:44:30Z")
Hi Andry I have faced a similar problem like you. But the elements of C matrix in my code are all functions of other field variables, like C(i,j) = Cij(m,n), here m and n are known field variables. And I use func C=[[C11,C12],[C21,C22] ]; to define C, but the same error of yours appears. Could you please lend me a hand with this? Thank you very much. Lee ---
**Author:** ![Lee](https://avatars.discourse-cdn.com/v4/letter/l/bc8723/32.png) [@Lee](https://community.freefem.org/u/Lee)\ **Post date:** [July 10, 2020, 3:57am UTC](https://community.freefem.org/t/matrix-vector-multiplication-in-weak-formulation/449/5 "2020-07-10T03:57:09Z")
Dear prj would you help me out with the similar problem? the short description is presented in this discussion part. And by the way, the varf equation I use is written in a very long formate and it takes a very long time to finish the computation of matrix A=equ(Th,Th). I want to retrive its matrix calculation form to speed up it. Could you also give me some suggestions on how to speed up the process? Thanks a lot. ---
**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\ **Post date:** [July 10, 2020, 6:12am UTC](https://community.freefem.org/t/matrix-vector-multiplication-in-weak-formulation/449/6 "2020-07-10T06:12:11Z")
You cannot use a `func`, use a `macro` instead. For speeding up your program, just use PETSc, see one of the many [examples](https://github.com/FreeFem/FreeFem-sources/blob/develop/examples/hpddm/README.md) + [this](http://jolivet.perso.enseeiht.fr/FreeFem-tutorial/#pf87) tutorial. Basically, replace ```auto mesh Th; func Pk = P1; varf vPb(u, v) = int2d(Th)(...); fespace Vh(Th, Pk); matrix A = equ(Vh, Vh); real[int] rhs = equ(0, Vh); real[int] prod = A * rhs; set(A, solver = sparsesolver); real[int] sol = A^-1 * rhs; ``` By ```auto mesh Th; func Pk = P1; load "PETSc" Mat A; createMat(Th, A, Pk) varf vPb(u, v) = int2d(Th)(...); fespace Vh(Th, Pk); Mat A = equ(Vh, Vh); real[int] rhs = equ(0, Vh); real[int] prod = A * rhs; set(A, sparams = "-pc_type lu"); real[int] sol = A^-1 * rhs; ``` You should get near-linear speedup for the assembly phase, and super-linear speedup for the solution phase (at least up until ~4 to 8 processes). ---
**Author:** ![Lee](https://avatars.discourse-cdn.com/v4/letter/l/bc8723/32.png) [@Lee](https://community.freefem.org/u/Lee)\ **Post date:** [July 13, 2020, 9:58am UTC](https://community.freefem.org/t/matrix-vector-multiplication-in-weak-formulation/449/7 "2020-07-13T09:58:03Z")
Dear prj Thanks for your reply, I will try PETSc. ---
**Author:** ![Lee](https://avatars.discourse-cdn.com/v4/letter/l/bc8723/32.png) [@Lee](https://community.freefem.org/u/Lee)\ **Post date:** [July 15, 2020, 2:57pm UTC](https://community.freefem.org/t/matrix-vector-multiplication-in-weak-formulation/449/8 "2020-07-15T14:57:54Z")
Dear prj My problem is a constrained non-linear problem, It works well using IPOPT provided in FreeFem++. But I have not found the parallel version of IPOPT in FreeFem++ so it take a very long time in sequence form, especially the formation of matrix HJ . Now I have exactly the form of J, DJ and HJ, can my problem be able to transformed using PETSc? Can you give me a general guideline for this transform process? Thanks a lot. ---
**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\ **Post date:** [July 15, 2020, 3:03pm UTC](https://community.freefem.org/t/matrix-vector-multiplication-in-weak-formulation/449/9 "2020-07-15T15:03:37Z")
Sure! You can use Tao, the Toolkit for Advanced Optimization, which is part of PETSc. The interface ressembles what you’d do with Ipopt, but of course, everything is parallel 🙂 You can have a look at one of these examples [[1](https://github.com/FreeFem/FreeFem-sources/blob/develop/examples/hpddm/minimal-surface-Tao-2d-PETSc.edp),[2](https://github.com/FreeFem/FreeFem-sources/blob/develop/examples/hpddm/orego-Tao-PETSc.edp),[3](https://github.com/FreeFem/FreeFem-sources/blob/develop/examples/hpddm/toy-Tao-PETSc.edp)]. Also, `example12.edp` from [this](http://jolivet.perso.enseeiht.fr/FreeFem-tutorial/main.pdf#page=174) tutorial. ---
**Author:** ![Lee](https://avatars.discourse-cdn.com/v4/letter/l/bc8723/32.png) [@Lee](https://community.freefem.org/u/Lee)\ **Post date:** [July 17, 2020, 10:50am UTC](https://community.freefem.org/t/matrix-vector-multiplication-in-weak-formulation/449/10 "2020-07-17T10:50:02Z")
Thank you very much. 😅 ---
**Author:** ![Lee](https://avatars.discourse-cdn.com/v4/letter/l/bc8723/32.png) [@Lee](https://community.freefem.org/u/Lee)\ **Post date:** [July 21, 2020, 9:42am UTC](https://community.freefem.org/t/matrix-vector-multiplication-in-weak-formulation/449/11 "2020-07-21T09:42:26Z")
Dear prj I have tried to read the example [https://github.com/FreeFem/FreeFem-sources/blob/develop/examples/hpddm/minimal-surface-Tao-2d-PETSc.edp](https://github.com/FreeFem/FreeFem-sources/blob/develop/examples/hpddm/minimal-surface-Tao-2d-PETSc.edp) using TAO and run it on my PC. I am confused by the problems below: 1. The output massege always prints “Solver terminated: -6 Line Search Failure” after several iterations. Does it mean the solving process is unsuccessful? If so, what should I do to make the line search successful? 2. Where can I found the detail explaination of the meaning of some functions like “changeNumbering”? I am not sure my guess is ture or not. This makes me confused on retrieving the desired global field variables. 3. My own problem is a little bit different from normal constrained nonlinear optimization, DJ is not the direct partial derivative of J and can not form a exact corresponding J. In my original IPOPT code, I must turn off the linesearch function (with the code “linesearch=false” in IPOPT) to get the results or it will trapped in an endless line search process. How should I deal with this problem with TAO? Thanks for your reply. ---
**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\ **Post date:** [July 21, 2020, 10:39am UTC](https://community.freefem.org/t/matrix-vector-multiplication-in-weak-formulation/449/12 "2020-07-21T10:39:21Z")
1. sorry, there was a bug in `minimal-surface-Tao-2d-PETSc.edp`, I fixed it in [this](https://github.com/FreeFem/FreeFem-sources/commit/b0caea3dea59b4a5061e0358f334ed92bc744f71) commit. With quasi-Newton, you should now get `Solver terminated: -2 Maximum Iterations`, and with Newton `Solution converged: ||g(X)|| <= gatol`. 2. please have a look a section 8 of [this](http://jolivet.perso.enseeiht.fr/FreeFem-tutorial/) tutorial. Also, you can start with a sequential example, and we’ll parallelize it afterwards. There is no need to bother with more than one process in the beginning. 3. you can try one of [these](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Tao/TaoType.html) methods, surely all of them don’t do linesearch. If you formulate exactly your problem, I may gave you a better answer. ---
**Author:** ![Lee](https://avatars.discourse-cdn.com/v4/letter/l/bc8723/32.png) [@Lee](https://community.freefem.org/u/Lee)\ **Post date:** [July 24, 2020, 7:17am UTC](https://community.freefem.org/t/matrix-vector-multiplication-in-weak-formulation/449/13 "2020-07-24T07:17:17Z")
Dear prj Thank you for your timely reply. I have managed to greatly simplify my problem and I tried to rewrite it in the same form of [[https://github.com/FreeFem/FreeFem-sources/blob/develop/examples/hpddm/minimal-surface-Tao-2d-PETSc.edp](https://github.com/FreeFem/FreeFem-sources/blob/develop/examples/hpddm/minimal-surface-Tao-2d-PETSc.edp)]. The original problem using IPOPT and the rewrite one are uploaded as follows. [Original\_Ipopt.edp](https://community.freefem.org/uploads/short-url/75Yh9E8Gu4TcHSYQImxcAuzbA8g.edp) (3.1 KB) [TAO\_test.edp](https://community.freefem.org/uploads/short-url/4vcVi7gxm3U3e0X6PkcPuV4empV.edp) (4.3 KB) The problem I want to solve is provided in Original\_Ipopt.edp. It is a simplified bound constrained-nonlinear-time dependent-multiple variable problem. It works fine using IPOPT. But I have not found a TAO example that meet all the requirement I need. So I tried to rewrite it using TAO as in TAO\_test.edp. In the TAO\_test.edp I uploaded, an error is encountered in the line 28 as createMat(Th, H, Pk); The message says that “meshN The Identifier meshN does not exist”. How to correct it? Could you please point out the incorrect parts in the TAO\_test.edp and give me your valuable suggestions on how to correct them? Expecting your reply. ---
**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\ **Post date:** [July 24, 2020, 8:35am UTC](https://community.freefem.org/t/matrix-vector-multiplication-in-weak-formulation/449/14 "2020-07-24T08:35:18Z")
Here, I’ve fixed part of the script: 1. macro `def` and `init` have to be scoped 2. macro `dimension` must be defined before including `macro_ddm.idp` 3. can’t use both `plotMPI` and `fespace Xh` (I’ll fix that in FreeFEM) 4. some functions were not defined at the right place The code runs, but I don’t know if it’s working properly. [TAO\_test.edp](https://community.freefem.org/uploads/short-url/6DR9SkigjxIgiNfbwWDqnSqdwXm.edp) (4.3 KB) ---
**Author:** ![Lee](https://avatars.discourse-cdn.com/v4/letter/l/bc8723/32.png) [@Lee](https://community.freefem.org/u/Lee)\ **Post date:** [July 24, 2020, 9:02am UTC](https://community.freefem.org/t/matrix-vector-multiplication-in-weak-formulation/449/15 "2020-07-24T09:02:45Z")
Thank you very much 😃 ---
**Author:** ![Lee](https://avatars.discourse-cdn.com/v4/letter/l/bc8723/32.png) [@Lee](https://community.freefem.org/u/Lee)\ **Post date:** [July 24, 2020, 10:49am UTC](https://community.freefem.org/t/matrix-vector-multiplication-in-weak-formulation/449/16 "2020-07-24T10:49:22Z")
Dear prj I have run the TAO\_test.edp, but the it stucked at line 36 and 37 as Line 36 plot(part, fill=1, value=1,wait=1,cmm=“part befor creatPartition”); Line 37 createPartition(Th, part[], P0) These two are used in the build up block of ThNo. The output is as below: — partition of unity built (in 2.956960e-02) — global numbering created (in 3.081115e-03) — global CSR created (in 4.659305e-03) Warning: May be a bug in your script, a part of the plot is wrong t (mesh or FE function, curve) =\> skip the item 1 in plot command Error plot item empty 0 Error plot item empty 0 Error plot item empty 0 Error plot item empty 0 Error plot item empty 0 Error plot item empty 0 Error plot item empty 0 Error plot item empty 0 Error plot item empty 0 Error plot item empty 0 Error plot item empty 0 Error plot item empty 0 Error plot item empty 0 Error plot item empty 0 Error plot item empty 0 The message “Error plot item empty 0” keep showing up and the next plot function in line 37did not happen. Could you please give me hint on how to fix it? Thank you very much. ---
**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\ **Post date:** [July 24, 2020, 11:28am UTC](https://community.freefem.org/t/matrix-vector-multiplication-in-weak-formulation/449/17 "2020-07-24T11:28:58Z")
Sorry, you need to add `buildDmesh(Th)` after the plot line 21. Then, `Th` will be distributed and both `createMat` and `createPartition` will work as expected. ---
**Author:** ![Lee](https://avatars.discourse-cdn.com/v4/letter/l/bc8723/32.png) [@Lee](https://community.freefem.org/u/Lee)\ **Post date:** [July 24, 2020, 11:36am UTC](https://community.freefem.org/t/matrix-vector-multiplication-in-weak-formulation/449/18 "2020-07-24T11:36:50Z")
It works! And the results seem to be the same as the Ipopt one I uploaded and the speed is extremely fast! Thank you very much! Do you have a full tutorial for the TAO in FreeFem++ to explain the functions like “createPartition”? The “Miscellaneous” part is FreeFem++ Doc is not clear enough for me when compared with explaination of IPOPT. Thank you again. ---
**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\ **Post date:** [July 24, 2020, 11:48am UTC](https://community.freefem.org/t/matrix-vector-multiplication-in-weak-formulation/449/19 "2020-07-24T11:48:34Z")
Tutorial about FreeFEM + parallelism is available [there](http://jolivet.perso.enseeiht.fr/FreeFem-tutorial/), section 8, read the PDF about PETSc, there are some explanation. If something is not clear, please let me know. List of Tao solvers are available [there](https://www.mcs.anl.gov/petsc/petsc-current/docs/tao_manual.pdf#page=35), if you see something you like that is not interfaced yet in FreeFEM, please let me know. ---
**Author:** ![Lee](https://avatars.discourse-cdn.com/v4/letter/l/bc8723/32.png) [@Lee](https://community.freefem.org/u/Lee)\ **Post date:** [July 26, 2020, 11:57am UTC](https://community.freefem.org/t/matrix-vector-multiplication-in-weak-formulation/449/20 "2020-07-26T11:57:09Z")
Dear prj I have implemented my complete problem using TAO, but the due to I can not find a correct optimizing function J for my DJ and HJ, almost each Newton iteration needs very long CG iteration computing J and DJ and can not reach convergence. In my original IPOPT program, I turned off the linesearch function and the J and DJ would be computed only once in each Newton iteration. In FreeFem++4.6 doc, it says that when linesearch is turned off, the method becomes a standard Newton algorithm instead of a primal-dual system(as in page 234). Would you please give me some suggestions on how to achieve this using TAO? I have tried all the possible TAO solvers for bound-constrained nonlinear problem and the problem continue to show up. I also tried the SNESSolve in PETSc, since it do not need a J. But I have questions on the meaning of parameters bPETSc and xPETSc in the SNESSolve implement line below: _SNESSolve(A, funcJ, funcRes, bPETSc, xPETSc,xl = xlPETSc, xu = xuPETSc, sparams = “-snes\_monitor -ksp\_converged\_reason -snes\_view -snes\_vi\_monitor -snes\_type vinewtonrsls -snes\_rtol 1.0e-6 -pc\_type lu”);_ It seems bPETSc is the initial value and xPETSc is the output, but I found that before entering SNESSolve, some use xPETSc=bPETSc, some use xPETSc=0 and some use xPETSc=0.1 in the provided examples using SNESSolve. Even for the same example, xPETSc=0.1 and xPETSc=0.2 would leads to different solutions. So would you please also tell me the exact meaning of this two parameters and their values? And by the way, the boundary condition is applied using the same way as TAO in a bound-constrained form, is this correct in SNESSolve or must be applied in the varf form? Thank you very much for your reply. [Next page](https://community.freefem.org/t/matrix-vector-multiplication-in-weak-formulation/449.md?page=2)