如何用C语言来实现快速PCA呢

由于需要求解的矩阵都是维度上万维，样本几十到几百个，之前用雅可比迭代法进行特征求解发i发现求解不出特征值和特征向量，如何用C语言来实现快速PCA呢

{
    int i, j, p, q, l;
    double fm, cn, sn, omega, x, y, d;
    l = 1;
    for (i = 0; i <= columns - 1; i++)    //赋值v=I 
    {
        v[i][i] = 1.0;
        for (j = 0; j <= columns - 1; j++)
            if (i != j) v[i][j] = 0;
    }
    while (true)
    {
        fm = 0.0;
        for (i = 0; i <= columns - 1; i++)
            for (j = 0; j <= i - 1; j++)
            {
            d = fabs(a[i][j]);
            if ((i != j) && (d>fm))   //取最大的a[p][q] 
            {
                fm = d; p = i; q = j;
            }
            }
        if (fm<eps) //如果达到给定的精度跳出循环 
        {
            return(1);
            break;
        }
        if (l>jt) //如果不能求解 
        {
            return(-1);
            break;
        }
        l = l + 1;     //记录迭代次数         
        x = -a[p][q];
        y = (a[q][q] - a[p][p]) / 2.0;
        omega = x / sqrt(x*x + y*y);
        if (y<0.0) omega = -omega;
        sn = 1.0 + sqrt(1.0 - omega*omega);
        sn = omega / sqrt(2.0*sn);  //正弦 
        cn = sqrt(1.0 - sn*sn);  //余弦 
        fm = a[p][p];
        a[p][p] = fm*cn*cn + a[q][q] * sn*sn + a[p][q] * omega;
        a[q][q] = fm*sn*sn + a[q][q] * cn*cn - a[p][q] * omega;
        a[p][q] = y*omega + a[p][q] * (2 * cn*cn - 1);
        a[q][p] = a[p][q];
        for (j = 0; j <= columns - 1; j++)
            if ((j != p) && (j != q))
            {
            fm = a[p][j];
            a[p][j] = fm*cn + a[q][j] * sn;
            a[q][j] = -fm*sn + a[q][j] * cn;
            }
        for (i = 0; i <= columns - 1; i++)
            if ((i != p) && (i != q))
            {
            fm = a[i][p];
            a[i][p] = fm*cn + a[i][q] * sn;
            a[i][q] = -fm*sn + a[i][q] * cn;
            }
        for (i = 0; i <= columns - 1; i++)
        {
            fm = v[i][p];
            v[i][p] = fm*cn + v[i][q] * sn;
            v[i][q] = -fm*sn + v[i][q] * cn;
        }
    }
}