Below Jacobi algorithm in povray format. I hope, it would be useful to somebody. As the example a ring is being drawn in the least-square plane defined by 5 points Pts[0..4]. The coordinates of points need to be initialized in the 3rd line. Of course, if the count is other than 5, you need to change variable NoPts in 2nd line.
#declare ooo=<0,0,0>;
#declare NoPts=5;
#declare Pts=array[NoPts]{p1,p2,p3,p4,p5}
// compute centroid
#declare Centroid=ooo;
#declare i = 0;
#while(i < NoPts )
#declare Centroid = Centroid+Pts[i];
#declare i = i + 1;
#end
#declare Centroid = Centroid / NoPts;
// Move points to centroid
#declare nPts=array[NoPts]{ooo,ooo,ooo,ooo,ooo}
#declare i = 0;
#while(i < NoPts )
#declare nPts[i] = Pts[i] - Centroid;
#declare i = i + 1;
#end
#declare SM=array[3][3]{ {0,0,0}, {0,0,0}, {0,0,0} } // 2nd momentum matrix for coordinates moved to the centroid
#declare EV=array[3][3]{ {1,0,0}, {0,1,0}, {0,0,1} } // eigen vectors
#declare i = 0;
#while(i < NoPts )
#declare SM[0][0] = nPts[i].x * nPts[i].x + SM[0][0];
#declare SM[0][1] = nPts[i].x * nPts[i].y + SM[0][1];
#declare SM[0][2] = nPts[i].x * nPts[i].z + SM[0][2];
#declare SM[1][1] = nPts[i].y * nPts[i].y + SM[1][1];
#declare SM[1][2] = nPts[i].y * nPts[i].z + SM[1][2];
#declare SM[2][2] = nPts[i].z * nPts[i].z + SM[2][2];
#declare i = i + 1;
#end
#declare SM[1][0] = SM[0][1];
#declare SM[2][0] = SM[0][2];
#declare SM[2][1] = SM[1][2];
// =============== JACOBI rotation
#declare tol = 0.00000001; // tolerance
#declare iterMax = 8;
#declare summ = 0; // control sum
#declare l = 0; // ----------- sanity check
#while(l < 3 )
#declare m = 0;
#while(m < 3 )
#declare summ = summ + abs(SM[l][m]);
#declare m = m + 1;
#end
#declare l = l + 1;
#end
#if (summ < 0.00005)
#debug concat("=============== STH WRONG \n")
#end
#declare j=0;
#while(j < iterMax)
#declare ssum = 0;
#declare amax = 0;
#declare row=1;
#while(row < 3) // ------- loop over rows
#declare col=0;
#while(col < row) // ------- loop over columns
#declare atmp = abs(SM[row][col]);
#if(atmp > amax)
#declare amax = atmp;
#end
#declare ssum = ssum + atmp;
#if(atmp < amax * 0.1)
#declare col = 5;
//#end
#else
#if(abs(ssum / summ) > tol ) // --- only for "huge" elements ;-)
#declare th = atan(2 * SM[row][col] / (SM[col][col] - SM[row][row])) / 2;
#declare sin_th = sin(th);
#declare cos_th = cos(th);
#declare k=0;
#while(k < 3)
#declare tmp2 = SM[k][col];
#declare SM[k][col] = cos_th * tmp2 + sin_th * SM[k][row];
#declare SM[k][row] = -sin_th * tmp2 + cos_th * SM[k][row];
#declare tmp2 = EV[k][col];
#declare EV[k][col] = cos_th * tmp2 + sin_th * EV[k][row];
#declare EV[k][row] = -sin_th * tmp2 + cos_th * EV[k][row];
#declare k = k + 1;
#end
#declare SM[col][col] = cos_th * SM[col][col] + sin_th * SM[row][col];
#declare SM[row][row] = -sin_th * SM[col][row] + cos_th * SM[row][row];
#declare SM[col][row] = 0;
#declare k=0;
#while(k < 3)
#declare SM[col][k] = SM[k][col];
#declare SM[row][k] = SM[k][row];
#declare k = k + 1;
#end
#end // end of loop for huge elements
#end
#declare col = col + 1;
#end // end col
#declare row = row + 1;
#end //end row
#declare j = j + 1; // ' ------- next iteration if not bigger than iterMAX
#end
// =============== end JACOBI
// find EV with the smallest eigenvalue
#declare EVmin = 10000;
#declare k=0;
#while (k < 3)
#if(SM[k][k] < EVmin)
#declare EVmin = SM[k][k];
#declare EVindex = k;
#end
#declare k=k+1;
#end
// TEST - draw the ring
#declare ringRadius=1;
#declare w = < EV[0][EVindex], EV[1][EVindex], EV[2][EVindex] >;
object{
torus { ringRadius*.75, Dash_Radius texture { Bond_Texture } }
Reorient_Trans(<0, 1, 0>,w)
translate Centroid
}