I think it's best to solve the individual 3 x 3 systems. The problems of doing such a large matrix, containing millions of 3 x 3 systems:
1. computational cost of assembly: assembling such a large matrix could be take a substantial amount of time (even with a sparse solver).
2. parallelization overheads: should there be a need to parallelize the code, the large matrix would have to be distributed over many processors and there would be communication overheads for passing information between processors. On the other hand, for individual 3 x 3 systems, the coefficients could easily be distributed over the various processors, after which it would be an "embarrassingly parallel" operation (i.e. no communication between processors would be necessary). 946
Furthermore, parallelization using OpenMP or MPI would not be too difficult. On the other hand, for the large matrix, it would probably be easier to make use of a library (such as PETSc: e.g. the example ex2.c is a nice start), rather than to write your own code.
Lastly, I would like to point out that for 3 x 3 systems, it's easy to work out an analytical formula using computer algebra systems (or using Wolfram Alpha). As an illustration, here's how it can be done with the Python package, sympy:
from sympy import *
x1,x2,x3=symbols('x1 x2 x3')
b1,b2,b3=symbols('b1 b2 b3')
a11,a12,a13=symbols('a11 a12 a13')
a21,a22,a23=symbols('a21 a22 a23')
a31,a32,a33=symbols('a31 a32 a33')
eq1=Eq(a11*x1+a12*x2+a13*x3,b1)
eq2=Eq(a21*x1+a22*x2+a23*x3,b2)
eq3=Eq(a31*x1+a32*x2+a33*x3,b3)
x=solve([eq1,eq2,eq3], [x1,x2,x3])
x.get(x1)
x.get(x2)
x.get(x3)
Just replace x1, x2, x3, b1, b2, b3, etc. with the actual variables in your existing code. If generation of C code is required:
[(c_name, c_code), (h_name, c_header)] = codegen(("x1",x.get(x1)), "C", "test", header=False, empty=False)
print(c_code)
[(c_name, c_code), (h_name, c_header)] = codegen(("x2",x.get(x2)), "C", "test", header=False, empty=False)
print(c_code)
[(c_name, c_code), (h_name, c_header)] = codegen(("x3",x.get(x3)), "C", "test", header=False, empty=False)
print(c_code)
Sample output (for just x1):
>>> print(c_code)
#include "test.h"
#include <math.h>
double x1(double a11, double a12, double a13, double a21, double a22, double a23, double a31, double a32, double a33, double b1, double b2, double b3) {
return (a12*a23*b3 - a12*a33*b2 - a13*a22*b3 + a13*a32*b2 + a22*a33*b1 - a23*a32*b1)/(a11*a22*a33 - a11*a23*a32 - a12*a21*a33 + a12*a23*a31 + a13*a21*a32 - a13*a22*a31);
}
And here it may make sense to pre-calculate the determinant to further reduce the computational cost, i.e.:
determinant = (a11*a22*a33 - a11*a23*a32 - a12*a21*a33 + a12*a23*a31 + a13*a21*a32 - a13*a22*a31); 904