If, as in the model problem, the grid and the domain are both rectangular, the operator produced by the FE method is similar to a finite difference operator. For this type of very regular problem it is probably wiser to use a finite difference method. However, the operator calculation code that I implemented does not assume this rectangular layout. The elements must be connected in a regular mesh, but they may have arbitrary shapes. This case is related to the conformally mapped finite difference scheme described by Fornberg (1989). Each method has some advantages. More time is spent in constructing the operator in the FE code, but boundary conditions are generally much easier to implement.
One common procedure for speeding up the CG method is to use a preconditioning operator. A diagonal preconditioner is obviously very easy to apply. If better preconditioners are desired the method depends on the layout of the data. If the global Galerkin operator is constructed, an incomplete Choleski preconditioner can be used. When the unassembled system consists of element matrices the preconditioning is usually done with element-by-element operators (Lee and Wathen, 1989). This is not simple to do when using the row-wise decomposition that I propose, but I believe that it should be possible to construct an efficient preconditioning operator with the row-wise data layout.