11 using namespace Eigen;
16 void iterJDI(
float64_t *C,
int *pMatSize,
int *pMatNumber,
int *ptn,
int *ptm,
20 double eps,
int itermax)
34 while (iter < itermax && ( (sh_max>eps) || (s_max>eps) ))
57 for (
int n=2;n<=*pMatSize;n++)
59 iterJDI(C, pMatSize, pMatNumber, &n, &m, s_max, sh_max, A);
65 for (
int i=0;i<MS;i++)
69 for (
int j=0;j<MS;j++)
70 col_norm[i] += A[i*MS+j]*A[i*MS+j];
72 col_norm[i] = sqrt(col_norm[i]);
78 for (
int i=0;i<MS;i++)
83 for (
int i=0;i<MS;i++)
84 d[i]=daux/col_norm[i];
86 for (
int i=0;i<MS;i++)
87 for (
int j=0;j<MS;j++)
90 for (
int k=0;k<MN;k++)
92 for (
int i=0;i<MS;i++)
94 for (
int j=0;j<MS;j++)
95 C[k*MS*MS+i*MS+j] *= (1/d[i])*(1/d[j]);
111 for (
int i = 0, d3 = 0; i < MN; i++, d3+=MS*MS)
113 tmm[i] = C[d3+m*MS+m];
114 tnn[i] = C[d3+n*MS+n];
115 tmn[i] = C[d3+n*MS+m];
121 for (
int i = 0; i < MN; i++)
123 G[i][0] = 0.5*(tmm[i]+tnn[i]);
124 G[i][1] = 0.5*(tmm[i]-tnn[i]);
128 for (
int i = 0; i < 3; i++)
130 for (
int j = 0; j <= i; j++)
134 for (
int k = 0; k < MN; k++)
135 GG[3*j+i] += G[k][i]*G[k][j];
137 GG[3*i+j] = GG[3*j+i];
141 for (
int i = 0; i < 3; i++)
146 EigenSolver<MatrixXd> eig;
149 MatrixXd eigenvectors = eig.pseudoEigenvectors();
150 VectorXd eigenvalues = eig.pseudoEigenvalueMatrix().diagonal();
152 eigenvectors.col(0) = eigenvectors.col(0).normalized();
153 eigenvectors.col(1) = eigenvectors.col(1).normalized();
154 eigenvectors.col(2) = eigenvectors.col(2).normalized();
156 memcpy(evectors, eigenvectors.data(), 9*
sizeof(
float64_t));
157 memcpy(evalues, eigenvalues.data(), 3*
sizeof(
float64_t));
160 if (fabs(evalues[1])<fabs(evalues[2]))
162 tmp_eval = evalues[1];
163 evalues[1] = evalues[2];
164 evalues[2] = tmp_eval;
165 memcpy(tmp_evec,&evectors[3],3*
sizeof(
float64_t));
166 memcpy(&evectors[3],&evectors[6],3*
sizeof(
float64_t));
167 memcpy(&evectors[6],tmp_evec,3*
sizeof(
float64_t));
169 if (fabs(evalues[0])<fabs(evalues[1]))
171 tmp_eval = evalues[0];
172 evalues[0] = evalues[1];
173 evalues[1] = tmp_eval;
174 memcpy(tmp_evec,evectors,3*
sizeof(
float64_t));
175 memcpy(evectors,&evectors[3],3*
sizeof(
float64_t));
176 memcpy(&evectors[3],tmp_evec,3*
sizeof(
float64_t));
178 if (fabs(evalues[1])<fabs(evalues[2]))
180 tmp_eval = evalues[1];
181 evalues[1] = evalues[2];
182 evalues[2] = tmp_eval;
183 memcpy(tmp_evec,&evectors[3],3*
sizeof(
float64_t));
184 memcpy(&evectors[3],&evectors[6],3*
sizeof(
float64_t));
185 memcpy(&evectors[6],tmp_evec,3*
sizeof(
float64_t));
191 for (
int i = 0; i < 9; i++)
192 aux[i] = evectors[i];
194 for (
int i = 0; i < 3; i++)
197 for (
int i = 0; i < 3; i++)
201 for (
int j = 0; j < 3; j++)
202 diagaux[i] += evectors[3*i+j] * aux[3*i+j];
204 diagaux[i] = 1/sqrt(fabs(diagaux[i]));
209 for (
int i = 0; i < 3; i++)
210 v[i] = evectors[3*ind_min+i]*diagaux[ind_min];
213 for (
int i = 0; i < 3; i++)
216 float64_t ch = sqrt( (1+sqrt(1+v[0]*v[0]))/2 );
218 float64_t c = sqrt( ( 1 + v[2]/sqrt(1+v[0]*v[0]) )/2 );
219 float64_t s = -v[1]/( 2*c*sqrt( (1+v[0]*v[0]) ) );
221 *sh_max =
std::max(*sh_max,fabs(sh));
231 for (
int i = 0; i < MS ;i++)
233 for (
int j = 0, d3 = 0; j < MN; j++, d3+=MS*MS)
235 h_slice1 = C[d3+i*MS+m];
236 h_slice2 = C[d3+i*MS+n];
237 buf[i][j] = rm11*h_slice1 + rm21*h_slice2;
238 C[d3+i*MS+n] = rm12*h_slice1 + rm22*h_slice2;
239 C[d3+i*MS+m] = buf[i][j];
243 for (
int i = 0; i < MS; i++)
245 for (
int j = 0, d3 = 0; j < MN; j++, d3+=MS*MS)
246 C[d3+m*MS+i] = buf[i][j];
248 for (
int i = 0; i < MS; i++)
250 for (
int j = 0, d3 = 0; j < MN; j++, d3+=MS*MS)
251 C[d3+n*MS+i] = C[d3+i*MS+n];
253 for (
int i = 0; i < MN; i++)
255 C[i*MS*MS+m*MS+m] = (rm11*rm11)*tmm[i]+(2*rm11*rm21)*tmn[i]+(rm21*rm21)*tnn[i];
256 C[i*MS*MS+n*MS+m] = (rm11*rm12)*tmm[i]+(rm11*rm22+rm12*rm21)*tmn[i]+(rm21*rm22)*tnn[i];
257 C[i*MS*MS+m*MS+n] = C[i*MS*MS+n*MS+m];
258 C[i*MS*MS+n*MS+n] = (rm12*rm12)*tmm[i]+(2*rm12*rm22)*tmn[i]+(rm22*rm22)*tnn[i];
264 for (
int i = 0; i < MS; i++)
268 A[m*MS+i] = rm22*col1 - rm12*col2;
269 A[n*MS+i] = -rm21*col1 + rm11*col2;
void sweepjedi(float64_t *C, int *pMatSize, int *pMatNumber, float64_t *s_max, float64_t *sh_max, float64_t *A)
all of classes and functions are contained in the shogun namespace
static SGMatrix< float64_t > diagonalize(SGNDArray< float64_t > C, SGMatrix< float64_t > V0=SGMatrix< float64_t >(NULL, 0, 0, false), double eps=CMath::MACHINE_EPSILON, int itermax=200)
Matrix::Scalar max(Matrix m)
void iterJDI(float64_t *C, int *pMatSize, int *pMatNumber, int *ptn, int *ptm, float64_t *s_max, float64_t *sh_max, float64_t *A)
static SGMatrix< T > create_identity_matrix(index_t size, T scale)