Dear Freefem community,
I have an issue with the addition of two matrices of type Mat. I read in the Forum that it should work since v4.9, however I receive the following error message :
Bmat = Amat + Tmat2; error operator + <N5PETSc14DistributedCSRIN5HPDDM7SchwarzINSt3__17complexIdEEEEEE>, <N5PETSc14DistributedCSRIN5HPDDM7SchwarzINSt3__17complexIdEEEEEE>
List of choices
Can you help me understand ?
Best,
Mayeul
Without the code, no.
Here is a relevant portion of my code, where I define the matrices I try to subtract, namely a variational formulation matrix and a spherical DtN matrix.
// ---------- DtN ----------
complex[int, int] SHvalues = SphericalHarmonics(maxDegree, xvals, yvals, zvals);
cout << “Spherical Harmonics computed” << endl;
int nbHarmonics = numberOfHarmonics(maxDegree);
int nbLegendre = numberOfLegendre(maxDegree);
complex[int,int] Earray(Vh.ndof, nbHarmonics);
complex[int] tmp(Vh.ndof);
matrix E;
VhSext sphericalHarmonic;
for (int l=0; l<=maxDegree; l++) {
sphericalHarmonic[] = SHvalues(:, numberLegendre(l,0));
cout << "Computing DtN coefficients for l = " << l << “, m = 0” << endl;
varf FourierCoef(u,v) = int2d(Th,10)(v*sphericalHarmonic);
Earray(:,numberHarmonic(l,0)) = FourierCoef(0,Vh);
for (int m=1; m<=l; m++) {
sphericalHarmonic\[\] = SHvalues(:, numberLegendre(l,m));
cout << "Computing DtN coefficients for l = " << l << ", m = +/-" << m << endl;
varf FourierCoef(u,v) = int2d(Th,10)(v*sphericalHarmonic);
tmp = FourierCoef(0,Vh);
Earray(:,numberHarmonic(l,-m)) = tmp;
Earray(:,numberHarmonic(l,m)) = conj(tmp);
}
}
E = Earray;
cout << “Matrix of spherical harmonics projections defined” << endl;
// Initialization of the diagonal matrix of the DtN coefficients
matrix D;
complex[int] Darray(nbHarmonics);
int n;
complex[int] valuesHankel = Hankel(maxDegree, k*Rext);
for (int l = 0; l <= maxDegree; l++) {
for (int m = -l; m <= l; m++) {
n = numberHarmonic(l,m);
Darray\[n\] = (l/Rext - k\*valuesHankel\[l + 1\]/valuesHankel\[l\])/(Rext\*Rext);
}
}
D = [Darray];
// cout << " D = " << D << endl;
cout << “Diagonal matrix of DtN coefficients defined” << endl;
// Final assembly
matrix T;
Mat Tmat;
T = E*D;
Tmat = T;
cout << “DtN operator defined” << endl;
Mat Eprime;
matrix Estar = E’;
Eprime = Estar;
MatMatMult(Tmat,Eprime,Tmat);
cout << “DtN operator computed” << endl;
// ---------- Verification of the DtN: the fundamental solution Phi(z), on a hard sphere ----------
matrix A;
Mat Amat, Bmat;
// ---------- Variationnal formulation ----------
varf a(u,v) =
int3d(Th)( dx(u) * dx(v)’ + dy(u) * dy(v)’ + dz(u) * dz(v)’ )
+ int3d(Th)(- k^2 * u * v’);
A = a(Vh, Vh);
Amat = A;
Bmat = Amat - Tmat;
cout << “Matrix assembled” << endl;
As I wrote above, the error I get is :
Bmat = Amat - Tmat; error operator - <N5PETSc14DistributedCSRIN5HPDDM7SchwarzINSt3__17complexIdEEEEEE>, <N5PETSc14DistributedCSRIN5HPDDM7SchwarzINSt3__17complexIdEEEEEE>
Thank you in advance for your help,
Mayeul
Your code is unreadable, make sure it’s inside a code block. On top of that, you can see in += and *= for PETSc Mat. · FreeFem/FreeFem-sources@63a93ac · GitHub that there is support for +=, not +. You should probably just do operations with FreeFEM matrix and then feed the result into a Mat.
(sorry I deleted my post by mistake)
Do you want to compile my code ? Because if so you might need the whole code with the definition of the functions that compute the Spherical harmonics and the spherical Hankel functions. If not, here is the code in a code box.
// ---------- DtN ----------
complex[int, int] SHvalues = SphericalHarmonics(maxDegree, xvals, yvals, zvals);
cout << “Spherical Harmonics computed” << endl;
int nbHarmonics = numberOfHarmonics(maxDegree);
int nbLegendre = numberOfLegendre(maxDegree);
complex[int,int] Earray(Vh.ndof, nbHarmonics);
complex[int] tmp(Vh.ndof);
matrix<complex> E;
VhSext sphericalHarmonic;
for (int l=0; l<=maxDegree; l++) {
sphericalHarmonic[] = SHvalues(:, numberLegendre(l,0));
cout << "Computing DtN coefficients for l = " << l << “, m = 0” << endl;
varf FourierCoef(u,v) = int2d(Th,10)(v*sphericalHarmonic);
Earray(:,numberHarmonic(l,0)) = FourierCoef(0,Vh);
for (int m=1; m<=l; m++) {
sphericalHarmonic\[\] = SHvalues(:, numberLegendre(l,m));
cout << "Computing DtN coefficients for l = " << l << ", m = +/-" << m << endl;
varf FourierCoef(u,v) = int2d(Th,10)(v*sphericalHarmonic);
tmp = FourierCoef(0,Vh);
Earray(:,numberHarmonic(l,-m)) = tmp;
Earray(:,numberHarmonic(l,m)) = conj(tmp);
}
}
E = Earray;
cout << “Matrix of spherical harmonics projections defined” << endl;
// Initialization of the diagonal matrix of the DtN coefficients
matrix<complex> D;
complex[int] Darray(nbHarmonics);
int n;
complex[int] valuesHankel = Hankel(maxDegree, k*Rext);
for (int l = 0; l <= maxDegree; l++) {
for (int m = -l; m <= l; m++) {
n = numberHarmonic(l,m);
Darray[n] = (l/Rext - k*valuesHankel[l + 1]/valuesHankel[l])/(Rext*Rext);
}
}
D = [Darray];
// cout << " D = " << D << endl;
cout << “Diagonal matrix of DtN coefficients defined” << endl;
// Final assembly
matrix<complex> T;
Mat<complex> Tmat;
T = E*D;
Tmat = T;
cout << “DtN operator defined” << endl;
Mat<complex> Eprime;
matrix<complex> Estar = E’;
Eprime = Estar;
MatMatMult(Tmat,Eprime,Tmat);
cout << “DtN operator computed” << endl;
// ---------- Verification of the DtN: the fundamental solution Phi(z), on a hard sphere ----------
matrix<complex> A;
Mat<complex> Amat, Bmat;
// ---------- Variationnal formulation ----------
varf a(u,v) =
int3d(Th)( dx(u) * dx(v)’ + dy(u) * dy(v)’ + dz(u) * dz(v)’ )
+ int3d(Th)(- k^2 * u * v’);
A = a(Vh, Vh);
Amat = A;
Bmat = Amat - Tmat;
Thank you for your help.
Sincerely,
Mayeul
Thank you very much. Indeed, I tried to do the operations with Freefem Matrix but the product MatMatMult(Tmat,Eprime,Tmat) is becoming very long when the meshsize is decreasing, and it is significantly shorter in Mat format. If you know a way to cast back the result of this product from Mat<complex> to matrix<complex> it could solve my issues.
I solved my issue with the += operator and the MatMatMult function thank you for your help.