A Slow Running Code

As a part of my project, I am using VPI's powermod function for computing matrix raised to a high power. I do know that powermod is faster than MATLAB's power function but my code runs very slowly (takes hours for a matrix of 512x512) and I cannot think of any way to fix it. How to use for-loop efficiently or how to use vectorization to make the code run faster. Any suggestion will be much appreciated. Here is the part of the code:

clc
clear

p = 321321197  
q = 321321211

p = vpi(p);
q = vpi(q);
N = p*q;
N = vpi(N);
N2 = vpi(N)^2;
g = N+1;
g = vpi(g);

I = magic(512);
[R, C] = size(I);

test = vpi(zeros(R, C));

for idx = 1 : R
    for idy = 1 : C
        test(idx, idy) = powermod(vpi(I(idx, idy)), N, N2);
    end
end

[KSSV]

First, you need to run a profiler and check where maximum time is spent. Read about profiling.

[Walter Roberson]

What happens if you skip the loops and try

test = powermod(vpi(I), N, N2)

[Mohsin Shah]

I think powermod does not support matrix input. I tried test = powermod(vpi(I), N, N2) and it gave this error: "Operands to the and && operators must be convertible to logical scalar values. Error in powermod (line 46) if (a > 2^26) (n > 2^26)"