# Matrix ^-1 \* matrix

**URL:** <https://community.freefem.org/t/matrix-1-matrix/1913>\
**Category:** General Discussion\
**Created:** [July 26, 2022, 2:52pm UTC](https://community.freefem.org/t/matrix-1-matrix/1913 "2022-07-26T14:52:01Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![murea](https://avatars.discourse-cdn.com/v4/letter/m/77aa72/32.png) [@murea](https://community.freefem.org/u/murea)\
**Post date:** [July 26, 2022, 2:52pm UTC](https://community.freefem.org/t/matrix-1-matrix/1913/1 "2022-07-26T14:52:01Z")

</div>

Hallo,

In order to compute the hessienne matrix for IPOPT, I want to compute  
BB1’ \* AA^-1 \* BB1  
where  
matrix AA=StokesForces3(Xh,Xh,solver=sparsesolver);  
matrix BB1=bb1(Mh,Xh);

If I try  
matrix HJ2;  
HJ2=AA^-1 \* BB1;  
HJ2 = BB1’ \* HJ2;  
the operation AA^-1 \* BB1 it is not accepted.

How, can I solve my problem, preferably, without using real[int,int] or lapack.

Thank you !  
Cornel

---

<div class="post-metadata">

**Author:** ![frederichecht](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/frederichecht/32/15_2.png) [@frederichecht](https://community.freefem.org/u/frederichecht)\
**Post date:** [July 27, 2022, 2:28pm UTC](https://community.freefem.org/t/matrix-1-matrix/1913/2 "2022-07-27T14:28:20Z")

</div>

Cornel,

the matrix BB1’ \* AA^-1 \* BB1 is full so no over possibility.

---

<div class="post-metadata">

**Author:** ![murea](https://avatars.discourse-cdn.com/v4/letter/m/77aa72/32.png) [@murea](https://community.freefem.org/u/murea)\
**Post date:** [July 27, 2022, 4:14pm UTC](https://community.freefem.org/t/matrix-1-matrix/1913/3 "2022-07-27T16:14:12Z")

</div>

Thanks Frédéric !

In fact AA=[[A, B’], [B, 0] ] and  
varf bb1([rh],[v1h,v2h,qh])=  
-int2d(Th)( (y1_v1h+y2_v2h)\*dHe(gtmp)\*rh/epsilon);

matrix BB1=bb1(Mh,Xh);

Xh [v1h,v2h,qh] are for Stokes  
Mh rh , gtmp is a kind of level set  
and dHe(gtmp) is non-zero only near in ring close the levelset=0

I have trayed with  
real[int,int] matBB1(Xh.ndof,Mh.ndof);  
matBB1 = 0.;  
int[int] II(1),JJ(1); real[int] CC(1);  
[II,JJ,CC]=BB1;

for(int i=0;i\<II.n;++i){  
matBB1(II(i),JJ(i)) = CC(i);  
};

real[int,int] matHJ2(Xh.ndof,Mh.ndof);

for(int j=0; j\<BB1.m; j++){  
matHJ2(:,j)=AA^-1 \* matBB1(:,j);  
}

matrix HJ2=matHJ2;  
HJ2 = BB1’ \* HJ2;  
HJ2 = 2\*HJ2;

but when I use real[int,int]  
I will limited to use coarse mesh ☹

Thanks a lot Frédéric !  
Cornel

---

<div class="post-metadata">

**Author:** ![frederichecht](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/frederichecht/32/15_2.png) [@frederichecht](https://community.freefem.org/u/frederichecht)\
**Post date:** [July 29, 2022, 11:57am UTC](https://community.freefem.org/t/matrix-1-matrix/1913/4 "2022-07-29T11:57:49Z")

</div>

> Remark, you can use a GMRES also to solve this linear system.

---

<div class="post-metadata">

**Author:** ![murea](https://avatars.discourse-cdn.com/v4/letter/m/77aa72/32.png) [@murea](https://community.freefem.org/u/murea)\
**Post date:** [July 29, 2022, 2:40pm UTC](https://community.freefem.org/t/matrix-1-matrix/1913/5 "2022-07-29T14:40:04Z")

</div>

Thanks Frédéric !

I have  
load “MUMPS\_seq”  
I put  
matrix AA=StokesForces3(Xh,Xh,solver=GMRES);  
in place of  
matrix AA=StokesForces3(Xh,Xh,solver=sparsesolver);  
but

for(int j=0; j\<BB1.m; j++){  
matHJ2(:,j)=AA^-1 \* matBB1(:,j);  
}  
takes a lot of time.

I put also  
matrix AA=StokesForces3(Xh,Xh,solver=LU,factorize=1);  
It is better, but I am forced to use in a unit square a mesh with h=1/40 when I use the HessianL  
IPOPT(J,dJ,HL,C,dC,gh,clb=clb,cub=cub,checkindex=1,structjacc=[gvi,gvj],maxiter=dk);

Without the HessianL (LBFGS)  
IPOPT(J,dJ,C,dC,gh,clb=clb,cub=cub,checkindex=1,structjacc=[gvi,gvj],maxiter=dk);  
I can use fine meshes h=1/200 etc  
For h=1/40, after 20 iterations, LBFGS obtains a smaller cost value than I use HessinL.  
Maybe, I have some errors when I compute dJ.

Thanks a lot Frédéric !  
Cornel