nomenclature (taken from memo of r. d. reasenberg dated
29 july 1975):
(v)   (c  f) (y)
rhs = u = ( ) = (    ) ( ) = bx,
(w)   (fa d) (z)
where fa = fadjoint
solutions are given by:
-1
y = (cbar)  *(vbar), where
-1                         -1
cbar = c - (f)*(d  )*(fa), vbar = v - (f)*(d  )*(w), and

storage arrays used in pprctl.  each is labeled with its use at the
successive operations: "<" = new assignment, "+" = modification,
"-" = use, "." = retention for later use, "=" = recovery of values,
"x" = scratch.

step:    s c f d i s p f  f  di ua da h v f  a e
e       n c t di di fa       d   di p n
t       v l r    w           r   fa r d
storage

ibuf6     < . . . . . .  .  .  .  .  . . -
ibuf7       < . . . . -  .  .  .  .  . . -
ibuf8             < . -
ibuf9   < . . . . . - <  .  -
ipartn  <
icrest  < . . . . . - .  .  .  .  .  -
idrest  < . . . . - - .  .  .  .  .  -
side      < . + - . - .  .  .  .  .  . -
vectk     < . + . . - +  .  .  .  .  -
b       x < < < + +   <  -  <  +  +  . . -  - -
* sav       x x x   x x                       x x

jptrc   <                                              1
jptrd   <                                              2
iuse            x                                      3
dvars             < . .  .  .  .  .  . . .  -          1-4

iptr    <                                              1
iptrc   < - -       =                                  2
iptrd   < . - -     = .  .  .  -  -                    3
scale           < -                                    1-4

* colfdi                x                                1-2
pvrow1          x              x  x  x                 1
vbar                                   <               1
f                                        <             1
pvrow2          x                    x                 2
u2ap                           <                       2
dap                               <                    2
c                                        <             2

mappings:
iptr(1-mparam) --> nparam
iptrc(1-mparam) --> ncparm
iptrd(1-mparam) --> ndparm
icrest(1-ncparm) --> nparam
idrest(1-ndparm) --> nparam
jptrc(1-nparam) --> ncparm
jptrd(1-nparam) --> ndparm

note that u (the c-part of rhs) can be restored using:
do 3000 i=1,ncparm
ia=iptr(icrest(i))
if(ia.le.0) go to 3000
side(ia)=side(ia)+sav(i)*weight
3000 continue