4
$\begingroup$

I’m going through some considerable effort to translate one of my codes from MATLAB. It’s a type of finite element code and I haven’t implemented the solver yet but comparing CPU times for simply assembling the global system between MATLAB and C, the latter is more than two orders of magnitude faster...

Anyways, I have no idea how to use PETSc but will attempt to learn it in the next few weeks. Just reading over some code and documentation, it seems like it might be pretty hard to translate my code for PETSc.

I noticed that there’s an option to read a matrix from a “file” (see here). I’m just wondering if I go through with translating my code as I’m doing now, can I just save matrix/rhs to a .dat file or something then read from there? I don’t have all that much experience with C so I don’t know if this would be difficult to do or not, especially for larger systems like I’ll have.

In my code, I’ll have to solve a linear system at every Newton iteration (it’s a nonlinear system), and combine this with the fact that the problem is time-dependent. The only time in my entire code that I’ll need PETSc is for the linear solve (admittedly, it’s a lot of linear solves...).

If I’m going to have to learn the syntax for PETSc now, I’d just like to know so I’m not frustrated later.

On another note, I’m really surprised how much faster the code is in C so far. Maybe my MATLAB coding isn’t that efficient, but I’ve been keeping the for loops at a minimum where possible.

$\endgroup$
6
$\begingroup$

I would recommend looking at a time-dependent example because PETSc can provide a lot more diagnostics if you formulate at that level. For example, you could use a Rosenbrock method, an additive Runge-Kutta IMEX method, or others. Some involve Newton iteration, but that is not required and is not always the best approach. To use those methods, you can start by just providing your residual function, $f(t,u)$, as in $\dot u = f(t,u)$, or alternatively, $g$ in $g(t,u,\dot u) = 0$. PETSc can automatically compute a Jacobian via coloring, and if you later decide to assemble the Jacobian for efficiency, PETSc can check whether what you assembled was correct and show you the errors. If you roll your own time integration and your own Newton, these algorithmic variants and debugging features will not be available unless you write them yourself.

The Users Manual explains the interfaces. Also see the manual pages for specific API documentation. The man pages are cross-referenced to examples. A good example to start with would be the Brusselator src/ts/examples/tutorials/ex25.c. You may also find slides from tutorials useful. If you get stuck, feel free to ask here or on the petsc-users mailing list.

Your proposal to use the file system will be a huge bottleneck even in serial. It complicates the workflow greatly and will limit the methods that you can use. But if you just want to test a matrix once, you can use PetscBinaryWrite() from MATLAB and load the system and solve using src/ksp/ksp/examples/tutorials/ex10.c.

$\endgroup$
  • $\begingroup$ Thanks, this was very helpful! I’m beginning to read up on all the details, but I’m just wondering in lay terms - if I assemble my system myself for now, and just want to use one of the linear solvers in PETSc (as opposed to LAPACK as I’m doing right now), will this be difficult given that my system is stored as an array of doubles right now? $\endgroup$ – Justin Dong Nov 20 '13 at 16:18
  • $\begingroup$ My understanding from your other question is that you have a sparse matrix. Storing that in a dense format turns an $O(n)$ algorithm into $O(n^2)$ storage, which is crippling if you want to solve larger problems. For testing, I recommend using PetscBinaryWrite() to write the sparse MATLAB matrix, plus a vector, to the file. You can use -ksp_view_solution binary:solution.petsc with the example (ex10.c), and read that with PetscBinaryRead() from MATLAB. $\endgroup$ – Jed Brown Nov 20 '13 at 16:49
  • $\begingroup$ ah, I see. I’m looking through some of the examples now. If I’m set on writing the assembly and such on my own for now, it seems that I could go through my code and translate the types for compatibility? i.e. int to PetscInt, double to PetscScalar, etc.? And same for implicit functions like malloc. I guess there’s nothing to it but to start coding a very small example. Also, I apologize for all these questions that must seem very silly; for someone just learning this gist of PETSc, the documentation is pretty dense. $\endgroup$ – Justin Dong Nov 20 '13 at 18:25
  • $\begingroup$ You don't have to use our typedefs, but if you solve problems larger than 2B degrees of freedom, you'll want to be able to switch to 64-bit integers, and if you have an extremely ill-conditioned system, you may want to use quad precision. Using PETSc's types will make those switches painless. You can use your own malloc(), but PETSc's does extra tracing in debug mode that is sometimes useful. $\endgroup$ – Jed Brown Nov 21 '13 at 6:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.