You will have to write a user macro to form L(neq,neqg) -> g(neq,neqg,1) = g(neq,neqg,2) where neqg is the number of RHS you want to solve and 'neq' the number of equations when you form TANG and is passed in cdata.h
Also A(neqg,neqg) must be formed by your user macro.
After you have formed the arrays call formhh(A,g,neqg,neq,neq)
The returned value in A will be what you want.
So to use you just use
TANG - which forms K and its factors
UMAC - where UMAC is replaced by the name of your user macro.
formhh uses the factors to do the multiple RHS solves -- it does do them 1 at a time.
I hope this is clear.