# Error during matrix-vector multiplication with ffddm **URL:** \ **Category:** General Discussion\ **Created:** [March 12, 2024, 1:25am UTC](https://community.freefem.org/t/error-during-matrix-vector-multiplication-with-ffddm/3030 "2024-03-12T01:25:58Z")\ **Posts on this page:** 3\ **Page:** 1
**Author:** ![luyu](https://avatars.discourse-cdn.com/v4/letter/l/c67d28/32.png) [@luyu](https://community.freefem.org/u/luyu)\ **Post date:** [March 12, 2024, 1:25am UTC](https://community.freefem.org/t/error-during-matrix-vector-multiplication-with-ffddm/3030/1 "2024-03-12T01:25:58Z")
Using ffddm, following the example `Maxwell-3d-simple.edp` in the `examples/ffddm` folder, I am trying to see what each operation does by printing out the results. So here I am multiplying the partition of unity matrix _pr#Dih[mpiRank(pr#commddm)]_ with a vector. I have defined x1 to be: ```auto complex[int] x1(rhs.n); x1 = 1; ``` The error lies in this line: ```auto x1 = MDih[mpiRank(Mcommddm)] * x1; ``` with the error: **Assertion fail : ((!this-\>N() || !this-\>M()) || !Ax.x.N() || !this-\>N() || &Ax.x[0] != &this-\>operator)** **line :560, in file ./…/femlib/RNM.hpp** **err code 6** I have checked that the output of MDih[mpiRank(Mcommddm)] .n and x1.n are the same so it shouldn’t be a problem of unmatching dimension. But I don’t really understand what the error is telling me. Below is my code for reference: ```auto load "iovtk" load "Element_Mixte3d" macro dimension 3// EOM include "ffddm.idp" load "msh3" macro def(i)[i, i#y, i#z]// EOM // vector field definition macro init(i)[i, i, i]// EOM // vector field initialization macro defPart(u)u// EOM // partition of unity definition macro initPart(u)u// EOM // partition of unity initialization macro Curl(ux, uy, uz)[dy(uz)-dz(uy), dz(ux)-dx(uz), dx(uy)-dy(ux)]// EOM macro CrossN(ux, uy, uz)[uy*N.z-uz*N.y, uz*N.x-ux*N.z, ux*N.y-uy*N.x]// EOM macro Curlabs(ux, uy, uz)[abs(dy(uz)-dz(uy)), abs(dz(ux)-dx(uz)), abs(dx(uy)-dy(ux))]// EOM func Pk = Edge03d; func PkPart = Edge03ds0; int Dirichlet = 1, Robin = 2; int ra = getARGV("-ra", 0); int ng = getARGV("-global", 30); int bc = getARGV("-bc", 1); real nk = getARGV("-k", 10.); real mul = getARGV("-mul", 40.); int[int] labs(1); if(bc == 1){labs = [Robin, Robin, Robin, Robin, Robin, Robin];} else{labs = [Dirichlet, Dirichlet, Dirichlet, Dirichlet, Dirichlet, Dirichlet];} // labs: list of surface labels, default=[1, 2, 3, 4, 5, 6] mesh3 Th = cube(ng, ng, ng, [x,y,z], label = labs); func f = exp(-((x-0.5)^2+(y-0.5)^2+(z-0.5)^2)*mul); int mysplit = 2; mesh3 Thc = cube(ng/mysplit, ng/mysplit, ng/mysplit, [x,y,z], label = labs); // mesh3 Thc = cube(7, 7, 7, [x,y,z], label = labs); // func k = 6*pi*nk; func k = nk; real epsilonprob = k^2;//0 k k^2; // epsilon of the problem real be = getARGV("-betaEprec",2.); // beta for epsilonEprec real epsilonEprec = k^be;//k^be; //macro Mdefplot(u)sqrt(real(u)^2+real(u#y)^2+real(u#z)^2)// vsym = 2; // symmetric, not hermitian vtgv = 1.e+30; vtgvelim = 1.e+30; ffddmbuildDmesh(M,Th,mpiCommWorld) ffddmbuildDfespaceEdge(M,M,complex,def,init,Pk,defPart,initPart,PkPart) macro Varf(varfName, meshName, PhName) varf varfName([Ex,Ey,Ez],[vx,vy,vz]) = int3d(meshName)(Curl(vx,vy,vz)'*Curl(Ex,Ey,Ez)) + int3d(meshName)(-(k^2-1i*epsilonprob)*[vx,vy,vz]'*[Ex,Ey,Ez]) + int2d(meshName,Robin)(1i*k*CrossN(vx,vy,vz)'*CrossN(Ex,Ey,Ez)) + on(Dirichlet,Ex=0,Ey=0,Ez=0); // EOM // for the preconditioner macro VarfEprec(varfName, meshName, PhName) varf varfName([Ex,Ey,Ez],[vx,vy,vz]) = int3d(meshName)(Curl(vx,vy,vz)'*Curl(Ex,Ey,Ez)) + int3d(meshName)(-(k^2-1i*epsilonEprec)*[vx,vy,vz]'*[Ex,Ey,Ez]) + int2d(meshName,Robin)(1i*k*CrossN(vx,vy,vz)'*CrossN(Ex,Ey,Ez)) + on(Dirichlet,Ex=0,Ey=0,Ez=0); // EOM // for the preconditioner macro VarfOpt(varfName, meshName, PhName) varf varfName([Ex,Ey,Ez],[vx,vy,vz]) = int3d(meshName)(Curl(vx,vy,vz)'*Curl(Ex,Ey,Ez)) + int3d(meshName)(-(k^2-1i*epsilonEprec)*[vx,vy,vz]'*[Ex,Ey,Ez]) + int2d(meshName,Robin)(1i*k*CrossN(vx,vy,vz)'*CrossN(Ex,Ey,Ez)) + int2d(meshName,10)(1i*k*CrossN(vx,vy,vz)'*CrossN(Ex,Ey,Ez)) + on(Dirichlet,Ex=0,Ey=0,Ez=0); // EOM macro Varfrhs(varfName, meshName, PhName) varf varfName([Ex,Ey,Ez],[vx,vy,vz]) = -int3d(meshName)([vx,vy,vz]'*[0,0,f]) + on(Dirichlet,Ex=0,Ey=0,Ez=0); // EOM ffddmsetup(M,M,Varf,Varfprec) if (ra == 1) {ffddmsetup(M, M, Varf, VarfOpt)} complex[int] rhs(1); ffddmbuildrhs(M,Varfrhs,rhs) complex[int] x0(rhs.n); x0 = 0; MVhi def(u); u[] = MfGMRES(x0, rhs, 1.e-6, 800, "right"); complex[int] x1(rhs.n); x1 = 1; cout << "Process " << mpiRank(Mcommddm) << "'s Di has " << MDih[mpiRank(Mcommddm)].n << " dofs." << endl; x1 = MDih[mpiRank(Mcommddm)] * x1; ``` ---
**Author:** ![phtournier](https://avatars.discourse-cdn.com/v4/letter/p/3da27b/32.png) [@phtournier](https://community.freefem.org/u/phtournier)\ **Post date:** [March 12, 2024, 9:05am UTC](https://community.freefem.org/t/error-during-matrix-vector-multiplication-with-ffddm/3030/2 "2024-03-12T09:05:50Z")
Unfortunately you cannot use the same vector as input and output for such operations in FreeFEM. You need to define a second vector as output : ```auto complex[int] x1(rhs.n), x2(rhs.n); x1 = 1; x2 = MDih[mpiRank(Mcommddm)] * x1; ``` or directly ```auto complex[int] x1(rhs.n); x1 = 1; complex[int] x2 = MDih[mpiRank(Mcommddm)] * x1; ``` ---
**Author:** ![luyu](https://avatars.discourse-cdn.com/v4/letter/l/c67d28/32.png) [@luyu](https://community.freefem.org/u/luyu)\ **Post date:** [March 12, 2024, 9:29am UTC](https://community.freefem.org/t/error-during-matrix-vector-multiplication-with-ffddm/3030/3 "2024-03-12T09:29:19Z")
Thanks! It works now. By the way, I found that the elementwise vector multiplication `x1 = x1 .* MDk[mpiRank(Mcommddm)]; ` works too. It seems a bit mysterious to me why the matrix-vector product won’t work. But thanks for the reply! It solved my problem.