Power Spectral Density estimate (PSD)
First of all, what do you mean when you say "Taken the first half of the resulting matrix (n/2+1)"?

The FFT will result in a vector, but maybe you mean taking only the first half of the FFT samples because the others are the negative frequencies and hence redundant?

Anyway, the process you describe is actually correct, and there is no need to do anything else. In your example with the resulting matrix a, you aparently took an FFT of size 8 (Or 16 if you then choped the second half of it). But what you do is simply take the average of all the first numbers of these vectors, and that would be an estimate of the Power spectral Density at f = 1*fs/NFFT. 
Taking the average of all the second number will be an estimate of PSD at f = 2*fs/NFFT, etc.

In the example the PDS would be

PSD = 1/3*[(1+2+4) , (2+5+5) ,(3+8+6) , ... ];

When you are analizing in real time as you carry out the computation, if the signal is Wide sense stationary, then it should converge or hover around some value, but if it is not, it will deviate over time to adjust to the new coming statistics of the signal.

If i understood you correctly, that should be it.

Edit:  oops i forgot about the important scaling factor. The scaling factor should be introduced to preserve the parseval theorem between the PSD and its transform (The autocorrelation). If you windowed your signal first, i think you should apply the scaling factor 1/sum(W*fs/NFFT), where W is your window vector. Since you are merely dividing by the integral of the window function in the frequency domain This is what you add to the estimate of the frequency component  when you do the FFT by windowing first. Because the FFT alone already preserves the energy.

Edit 2:

I am even more confused with the 2-D thing. Let me try an ask you a few questions:

1 - What type of signal are you trying to process? I.e: What type of application is it? Is it a unidimensional signal, like Audio sampled at 44.1Khz? Or is it a multidimensional signal like an image or video?

2 - Where are you performing these calculations? Matlab, Octave, python, or some other?

3 - When you say 2D vector (maybe could be less confusing to call this a Matrix), What does the numbers inside mean? what are the columns and what are the rows? 

4 - Did you check the link to mathworks? i think you will find what you are trying to do, is well explained there.

Also, just to make sure you understand the backgroun well enough, you could refer to some classical books that have the topic very well explained:


Digital Signal Processing, Principles, Algorithms And Applications -Proakis &amp; Manolakis - 3rd Ed. Ch. 12


And the Father of all books in spectral Estimation:


Spectral Analysis of Signals - Stoica &amp; Moses


Edit 3:

Ahh, ok, now i get what you are trying to do. So, here's the thing. I gave a quick read to the sections on FFT and PSD in the books by Stoica and Proakis, and still can't quite understand why they propose a scaling factor of simply 1/N to the FFT, while the matlab function for PSD using periodogram method scales the FFT by 2/(Fs*N), being N the length of the FFt, in your case 256.

So i don't know a lot of C++ but i would say do the following:

1 - Take first frame X, which consists of 256 samples of audio signal.

2 - apply 256 FFT to it. you obtain Xfft.

3 - To each complex value Xfft(i), for i = 0...255, do xfft(i)*conj(xfft(i))
which is equal to abs(xfft(i))^2. this you already did,so up to here nothing new.

4 - Multiply the whole vector Xfft by 2/(Fs*N).

This will give you the Power spectral density for your frame X. Now take the next frame and do the same, And when you reach a number of frames that you would like to plot, you stop and plot your vectors.

now since you are taking frames every 128 samples, your estimates at each time step will be sort of crammed up together, so a trick that youwould use is the fact that a Speech/audio signal is approximately wide sense statinary (WSS) in an interval of 20-50 msec (See Speech Enhancement - Jacob Benesty).
So what you could do, is assuming you have an Fs = 44100, that gives you a period of 50msec of around 2205 samples.

With the dimensions of the frames you took you could easily average 15 or 16 overlaped frames and get a Smoother estimate of the PSD. But the choice is yours, experiment with it and have fun!