Matrix ^-1 * matrix

General Discussion
1

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

2

Cornel,

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

3

Thanks Frédéric !

In fact AA=[ [A, B’], [B, 0] ] and
varf bb1([rh],[v1h,v2h,qh])=
-int2d(Th)( (y1v1h+y2v2h)*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 :frowning:

Thanks a lot Frédéric !
Cornel

4

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

5

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