GPULIB continued

It is often necessary to calculate the covariance matrix of a multispectral image from the array of observation vectors (pixels). The program

pro test_gpuMatrix_Multiply
; initialize
gpuinit
for i=3,6 do begin
; array of 6-element observation vectors (rows)
A = randomu(s,6,10L^i)
gpuPutArr,A,A_gpu
; calculate covariance matrix C on the CPU
start = systime(2)
for j=0L,99 do C = $
Matrix_Multiply(A,A,/btranspose)
CPUtime = systime(2)-start
print, 'Length : ',10L^i
print, 'CPU time: ',CPUtime
; now on the GPU
start = systime(2)
for j=0L,999 do $
gpuMatrix_Multiply,A_gpu,A_gpu,C_gpu,/btranspose
GPUtime = systime(2)-start
print,'GPU time: ', GPUtime/10
print,'Speedup : ', 10*CPUtime/GPUtime
gpuFree,A_gpu
endfor
; check that the results are the same
gpuGetArr,C_gpu,C1
print,'Check:'
print, determ(C),determ(C1)
gpuFree,C_gpu
end

produces

% Compiled module: TEST_GPUMATRIX_MULTIPLY.
Length : 1000
CPU time: 0.030999899
GPU time: 0.0062000036
Speedup : 4.9999808
Length : 10000
CPU time: 0.25000000
GPU time: 0.0046999931
Speedup : 53.191567
Length : 100000
CPU time: 0.96900010
GPU time: 0.0046999931
Speedup : 206.17054
Length : 1000000
CPU time: 9.5620000
GPU time: 0.0032000065
Speedup : 2988.1190
Check:
6.33251e+030 6.33251e+030

So for instance the covariance matrix of a 1000x1000x6 multispectral image could be calculated 3000 times faster on the GPU!? My ENVI/IDL change detection extension iterates on the weighted covariance matrix of a bitemporal image. Sooo ...