Sorry, the right hand side should be A^t*b:
Dune::FieldMatrix<Scalar, numSurfaceReactions, numSurfaceSpecies> A; Dune::FieldMatrix<Scalar, numSurfaceSpecies, numSurfaceReactions> At; for(int speciesIdx=0; speciesIdx<numSurfaceSpecies; speciesIdx++) for(int reactionIdx=0; reactionIdx<numSurfaceReactions; reactionIdx++) { A[reactionIdx][speciesIdx] = ... At[speciesIdx][reactionIdx] = ... } auto AtA = A.leftmultiply(At); // right hand side Dune::FieldVector<Scalar, numSurfaceReactions> b; for(int reactionIdx=0; reactionIdx<numSurfaceReactions; reactionIdx++) b[reactionIdx] = … Dune::FieldVector<Scalar, numSurfaceSpecies> Atb; A.mv(b, Atb); // solution Dune::FieldVector<Scalar, numSurfaceSpecies> x; AtA.solve(x, Atb); ________________________________ Von: Dumux <dumux-boun...@listserv.uni-stuttgart.de> im Auftrag von Flemisch, Bernd, apl.Prof.Dr. <bernd.flemi...@iws.uni-stuttgart.de> Gesendet: Montag, 11. Dezember 2017 10:31:07 An: DuMuX User Forum Betreff: Re: [DuMuX] Dune::BCRSMatrix In that case (numSurfaceReactions > numSurfaceSpecies), you have to go for a least-squares solution min_{x in R^numSurfaceSpecies} || A*x - b || which can be achieved by solving (A^t*A)*x = b I have not seen a dedicated ISTL solver for this. My indexing was wrong in the first place. You should modify to Dune::FieldMatrix<Scalar, numSurfaceReactions, numSurfaceSpecies> A; Dune::FieldMatrix<Scalar, numSurfaceSpecies, numSurfaceReactions> At; for(int speciesIdx=0; speciesIdx<numSurfaceSpecies; speciesIdx++) for(int reactionIdx=0; reactionIdx<numSurfaceReactions; reactionIdx++) { A[reactionIdx][speciesIdx] = ... At[speciesIdx][reactionIdx] = ... } auto AtA = A.leftmultiply(At); // right hand side Dune::FieldVector<Scalar, numSurfaceReactions> b; for(int reactionIdx=0; reactionIdx<numSurfaceReactions; reactionIdx++) b[reactionIdx] = … // solution Dune::FieldVector<Scalar, numSurfaceSpecies> x; AtA.solve(x, b); ________________________________ Von: Dumux <dumux-boun...@listserv.uni-stuttgart.de> im Auftrag von georg.fut...@dlr.de <georg.fut...@dlr.de> Gesendet: Montag, 11. Dezember 2017 10:10:20 An: dumux@listserv.uni-stuttgart.de Betreff: Re: [DuMuX] Dune::BCRSMatrix Hi Bernd, Actually, the system might be overdetermined, so FieldMatrix solve() will not necessarily work. Do I have any options to solve such a system with ISTL? Kind regards Georg Von: Dumux [mailto:dumux-boun...@listserv.uni-stuttgart.de] Im Auftrag von Flemisch, Bernd, apl.Prof.Dr. Gesendet: Montag, 11. Dezember 2017 10:03 An: DuMuX User Forum Betreff: Re: [DuMuX] Dune::BCRSMatrix Hi Georg, for the system that you describe, there is no need to wrap things in BCRSMatrices and BlockVectors. You can use directly FieldMatrices and -Vectors: Dune::FieldMatrix<Scalar, numSurfaceReactions, numSurfaceSpecies> A; for(int speciesIdx=0; speciesIdx<numSurfaceSpecies; speciesIdx++) for(int reactionIdx=0; reactionIdx<numSurfaceReactions; reactionIdx++) A[speciesIdx][reactionIdx] = ... // right hand side Dune::FieldVector<Scalar, numSurfaceReactions> b; for(int reactionIdx=0; reactionIdx<numSurfaceReactions; reactionIdx++) b[reactionIdx] = … // solution Dune::FieldVector<Scalar, numSurfaceSpecies> x; A.solve(x, b); I hope that numSurfaceReactions == numSurfaceSpecies? You would only need the more complicated setting if your system size was not known at compile time. Kind regards Bernd ________________________________ Von: Dumux <dumux-boun...@listserv.uni-stuttgart.de<mailto:dumux-boun...@listserv.uni-stuttgart.de>> im Auftrag von georg.fut...@dlr.de<mailto:georg.fut...@dlr.de> <georg.fut...@dlr.de<mailto:georg.fut...@dlr.de>> Gesendet: Montag, 11. Dezember 2017 09:29:07 An: dumux@listserv.uni-stuttgart.de<mailto:dumux@listserv.uni-stuttgart.de> Betreff: [DuMuX] Dune::BCRSMatrix Hello Dumux, I would like to solve a linear system of equations (A*x = b) using the Dune ISTL. To use the Dumux implementation (seqsolverbackend.hh), Dune::FieldVectors need to be wrapped in Dune::Blockvectors and a Dune::FieldMatrix is wrapped in a Dune::BCRSMatrix. I am not familiar with these data types and the Dune documentation on the homepage is really confusing. Since the system is not too large, I guess I can solve everything in a single block? This is the code I use right now: Dune::FieldMatrix<Scalar, numSurfaceReactions, numSurfaceSpecies> A; for(int speciesIdx=0; speciesIdx<numSurfaceSpecies; speciesIdx++) for(int reactionIdx=0; reactionIdx<numSurfaceReactions; reactionIdx++) assign values… // right hand side Dune::BlockVector<Dune::FieldVector<Scalar, numSurfaceReactions> > b; b.resize(1); // does this mean there is only one block containing a FieldVector ? for(int reactionIdx=0; reactionIdx<numSurfaceReactions; reactionIdx++) b[0][reactionIdx] = … // solution Dune::BlockVector<Dune::FieldVector<Scalar, numSurfaceSpecies> > x; x.resize(1); typedef typename Dune::FieldMatrix<Scalar, numSurfaceReactions, numSurfaceSpecies> MatrixBlock; typedef typename Dune::BCRSMatrix<MatrixBlock> ISTLMatrix; ISTLMatrix A2(1, 1, ISTLMatrix::random); // does this mean there is only one block containing a FieldMatrix ? for(int speciesIdx=0; speciesIdx<numSurfaceSpecies; speciesIdx++) for(int reactionIdx=0; reactionIdx<numSurfaceReactions; reactionIdx++) A2[0][0][reactionIdx][speciesIdx] = … HOW TO ASSIGN VALUES??? I get a segmentation fault from the last line. Is there better information on BCRS matrices somewhere? Any hints on how to set up this matrix are greatly appreciated! Best regards Georg Futter —————————————————————————— German Aerospace Center (DLR) Institute of Engineering Thermodynamics | Computational Electrochemistry | Pfaffenwaldring 38-40 | 70569 Stuttgart Dipl.-Ing. Georg Futter | Ph.D. student Telefon 0711/6862-8135 | georg.fut...@dlr.de<mailto:georg.fut...@dlr.de> www.DLR.de<http://www.dlr.de/>
_______________________________________________ Dumux mailing list Dumux@listserv.uni-stuttgart.de https://listserv.uni-stuttgart.de/mailman/listinfo/dumux