00001
00002 #ifdef VCL_NEEDS_PRAGMA_INTERFACE
00003 #pragma implementation
00004 #endif
00005
00006
00007
00008 #include "vnl_sparse_symmetric_eigensystem.h"
00009 #include "vnl_sparse_lu.h"
00010 #include <vnl/vnl_vector_ref.h>
00011 #include <vcl_cassert.h>
00012 #include <vcl_cstring.h>
00013 #include <vcl_cstdlib.h>
00014 #include <vcl_iostream.h>
00015 #include <vcl_vector.h>
00016
00017 #include <vnl/algo/vnl_netlib.h>
00018
00019 static vnl_sparse_symmetric_eigensystem * current_system = 0;
00020
00021 #ifdef VCL_SUNPRO_CC
00022 # define FUNCTION extern "C"
00023 #else
00024 # define FUNCTION static
00025 #endif
00026
00027
00028
00029
00030
00031 FUNCTION
00032 void sse_op_callback(const long* n,
00033 const long* m,
00034 const double* p,
00035 double* q)
00036 {
00037 assert(current_system != 0);
00038
00039 current_system->CalculateProduct(*n,*m,p,q);
00040 }
00041
00042
00043
00044
00045
00046
00047 FUNCTION
00048 void sse_iovect_callback(const long* n,
00049 const long* m,
00050 double* q,
00051 const long* j,
00052 const long* k)
00053 {
00054 assert(current_system != 0);
00055
00056 if (*k==0)
00057 current_system->SaveVectors(*n,*m,q,*j-*m);
00058 else if (*k==1)
00059 current_system->RestoreVectors(*n,*m,q,*j-*m);
00060 }
00061
00062 vnl_sparse_symmetric_eigensystem::vnl_sparse_symmetric_eigensystem()
00063 : nvalues(0), vectors(0), values(0)
00064 {
00065 }
00066
00067 vnl_sparse_symmetric_eigensystem::~vnl_sparse_symmetric_eigensystem()
00068 {
00069 delete[] vectors; vectors = 0;
00070 delete[] values; values = 0;
00071 for (unsigned i=0; i<temp_store.size(); ++i)
00072 delete temp_store[i];
00073 temp_store.clear();
00074 }
00075
00076
00077
00078
00079
00080
00081
00082
00083 int vnl_sparse_symmetric_eigensystem::CalculateNPairs(vnl_sparse_matrix<double>& M,
00084 int n,
00085 bool smallest,
00086 long nfigures)
00087 {
00088 mat = &M;
00089
00090
00091 if (vectors) {
00092 delete[] vectors; vectors = 0;
00093 delete[] values; values = 0;
00094 }
00095 nvalues = 0;
00096
00097 current_system = this;
00098
00099 long dim = mat->columns();
00100 long nvals = (smallest)?-n:n;
00101 long nperm = 0;
00102 long nmval = n;
00103 long nmvec = dim;
00104 vcl_vector<double> temp_vals(n*4);
00105 vcl_vector<double> temp_vecs(n*dim);
00106
00107
00108 long nblock = (dim<60) ? dim/6 : 10;
00109
00110
00111 long maxop = dim*10;
00112
00113
00114 long maxj = maxop*nblock;
00115 long t1 = 6*nblock+1;
00116 if (maxj < t1) maxj = t1;
00117 if (maxj < 40) maxj = 40;
00118
00119
00120
00121 int work_size = dim*nblock;
00122 int t2 = maxj*(2*nblock+3) + 2*n + 6 + (2*nblock+2)*(nblock+1);
00123 if (work_size < t2) work_size = t2;
00124 work_size += 2*dim*nblock + maxj*(nblock + n + 2) + 2*nblock*nblock + 3*n;
00125 vcl_vector<double> work(work_size+10);
00126
00127
00128 for (int i=0; i<dim*nblock; ++i)
00129 work[i] = 0.0;
00130
00131 vcl_vector<long> ind(n);
00132
00133 long ierr = 0;
00134
00135 v3p_netlib_dnlaso_(sse_op_callback, sse_iovect_callback,
00136 &dim, &nvals, &nfigures, &nperm,
00137 &nmval, &temp_vals[0],
00138 &nmvec, &temp_vecs[0],
00139 &nblock,
00140 &maxop,
00141 &maxj,
00142 &work[0],
00143 &ind[0],
00144 &ierr);
00145 if (ierr > 0)
00146 {
00147 if (ierr & 0x1)
00148 vcl_cerr << "Error: vnl_sparse_symmetric_eigensystem: N < 6*NBLOCK\n";
00149 if (ierr & 0x2)
00150 vcl_cerr << "Error: vnl_sparse_symmetric_eigensystem: NFIG < 0\n";
00151 if (ierr & 0x4)
00152 vcl_cerr << "Error: vnl_sparse_symmetric_eigensystem: NMVEC < N\n";
00153 if (ierr & 0x8)
00154 vcl_cerr << "Error: vnl_sparse_symmetric_eigensystem: NPERM < 0\n";
00155 if (ierr & 0x10)
00156 vcl_cerr << "Error: vnl_sparse_symmetric_eigensystem: MAXJ < 6*NBLOCK\n";
00157 if (ierr & 0x20)
00158 vcl_cerr << "Error: vnl_sparse_symmetric_eigensystem: NVAL < max(1,NPERM)\n";
00159 if (ierr & 0x40)
00160 vcl_cerr << "Error: vnl_sparse_symmetric_eigensystem: NVAL > NMVAL\n";
00161 if (ierr & 0x80)
00162 vcl_cerr << "Error: vnl_sparse_symmetric_eigensystem: NVAL > MAXOP\n";
00163 if (ierr & 0x100)
00164 vcl_cerr << "Error: vnl_sparse_symmetric_eigensystem: NVAL > MAXJ/2\n";
00165 if (ierr & 0x200)
00166 vcl_cerr << "Error: vnl_sparse_symmetric_eigensystem: NBLOCK < 1\n";
00167 }
00168 else if (ierr < 0)
00169 {
00170 if (ierr == -1)
00171 vcl_cerr << "Error: vnl_sparse_symmetric_eigensystem:\n"
00172 << " poor initial vectors chosen\n";
00173 else if (ierr == -2)
00174 vcl_cerr << "Error: vnl_sparse_symmetric_eigensystem:\n"
00175 << " reached maximum operations " << maxop
00176 << " without finding all eigenvalues,\n"
00177 << " found " << nperm << " eigenvalues\n";
00178 else if (ierr == -8)
00179 vcl_cerr << "Error: vnl_sparse_symmetric_eigensystem:\n"
00180 << " disastrous loss of orthogonality - internal error\n";
00181 }
00182
00183
00184 nvalues = n;
00185 vectors = new vnl_vector<double>[n];
00186 values = new double[n];
00187 for (int i=0; i<n; ++i) {
00188 values[i] = temp_vals[i];
00189 #if 0
00190 vcl_cout << "value " << temp_vals[i]
00191 << " accuracy " << temp_vals[i+n*2] << vcl_endl;
00192 #endif
00193 vnl_vector<double> vec(dim,0.0);
00194 for (int j=0; j<dim; ++j)
00195 vec[j] = temp_vecs[j + dim*i];
00196 vectors[i] = vec;
00197 }
00198
00199
00200 for (unsigned i=0; i<temp_store.size(); ++i)
00201 delete [] temp_store[i];
00202 temp_store.clear();
00203
00204 return ierr;
00205 }
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217 int vnl_sparse_symmetric_eigensystem::CalculateNPairs(
00218 vnl_sparse_matrix<double>& A, vnl_sparse_matrix<double>& B, int nEV,
00219 double tolerance, int numberLanczosVecs,
00220 bool smallest, bool magnitude,
00221 int maxIterations,
00222 double sigma)
00223 {
00224 mat = &A;
00225 Bmat = &B;
00226
00227
00228 if (vectors) {
00229 delete[] vectors; vectors = NULL;
00230 delete[] values; values = NULL;
00231 }
00232 nvalues = 0;
00233
00234 const long whichLength = 2;
00235 char which[whichLength + 1];
00236 which[whichLength] = '\0';
00237 if (smallest)
00238 which[0] = 'S';
00239 else
00240 which[0] = 'L';
00241
00242 if (magnitude)
00243 which[1] = 'M';
00244 else
00245 which[1] = 'A';
00246
00247 long matSize = mat->columns();
00248 long ido = 0;
00249
00250 long nconv = 0L;
00251 long numberLanczosVecsL = numberLanczosVecs;
00252 long nEVL = nEV;
00253
00254 double *resid = new double[matSize];
00255 vcl_memset((void*) resid, 0, sizeof(double)*matSize);
00256
00257 if (maxIterations <= 0)
00258 maxIterations = nEVL * 100;
00259
00260 if (numberLanczosVecsL <= 0)
00261 numberLanczosVecsL = 2 * nEVL + 1;
00262 numberLanczosVecsL = (numberLanczosVecsL > matSize ? matSize : numberLanczosVecsL);
00263 double *V = new double[matSize * numberLanczosVecsL + 1];
00264
00265 #define DONE 99
00266 const int genEigProblemLength = 1;
00267 char genEigProblem = 'G';
00268 long info = 0;
00269
00270 #define IPARAMSIZE 12
00271 long iParam[IPARAMSIZE];
00272
00273
00274 iParam[0] = 0;
00275 iParam[1] = 1;
00276 iParam[2] = 0;
00277 iParam[3] = maxIterations;
00278 iParam[4] = 1;
00279
00280
00281 iParam[5] = 0;
00282 iParam[6] = 0;
00283
00284 long mode;
00285
00286
00287
00288
00289 vnl_sparse_matrix<double> OP;
00290 if (sigma != 0.0)
00291 {
00292
00293
00294
00295 mode = 3;
00296
00297
00298 OP = B;
00299 OP *= sigma;
00300 OP = A - OP;
00301
00302 }
00303 else
00304 {
00305
00306
00307 mode = 2;
00308 OP = B;
00309 }
00310
00311 iParam[7] = mode;
00312
00313
00314 vnl_sparse_lu opLU(OP);
00315
00316
00317
00318 iParam[8] = 0;
00319
00320
00321 iParam[9] = 0;
00322 iParam[10] = 0;
00323 iParam[11] = 0;
00324
00325
00326
00327
00328 long iPntr[IPARAMSIZE];
00329 for (int clrIx = 0; clrIx < IPARAMSIZE; clrIx++)
00330 iPntr[clrIx]= 0;
00331
00332
00333 double *workd = new double[3 * matSize + 1];
00334
00335
00336 long lworkl = numberLanczosVecsL * (numberLanczosVecsL+9);
00337
00338
00339 double *workl = new double[lworkl + 1];
00340
00341
00342 bool basisCalculated = false;
00343
00344 vnl_vector<double> workVector;
00345
00346 while (!basisCalculated)
00347 {
00348
00349 v3p_netlib_dsaupd_(
00350 &ido, &genEigProblem, &matSize, which,
00351 &nEVL, &tolerance, resid, &numberLanczosVecsL, &V[1], &matSize,
00352 &iParam[1], &iPntr[1], &workd[1], &workl[1], &lworkl, &info,
00353 genEigProblemLength, whichLength);
00354
00355
00356 if (ido==DONE)
00357 {
00358 nconv = iParam[5];
00359 basisCalculated = true;
00360 break;
00361 }
00362 else
00363 {
00364 switch (info) {
00365 case -8:
00366 case -9:
00367 case -9999:
00368 return info;
00369 break;
00370 case 0:
00371 case 1:
00372 case 3:
00373 break;
00374 default :
00375 return info;
00376 }
00377
00378
00379 if (ido == -1)
00380 iPntr[3] = iPntr[2] + matSize;
00381
00382 vnl_vector_ref<double> x(matSize, &workd[iPntr[1]]);
00383 vnl_vector_ref<double> y(matSize, &workd[iPntr[2]]);
00384 vnl_vector_ref<double> z(matSize, &workd[iPntr[3]]);
00385
00386 switch (ido)
00387 {
00388 case -1:
00389
00390 if (mode != 2)
00391 B.mult(x, z);
00392
00393 case 1:
00394
00395 if (mode != 2)
00396 opLU.solve(z, &y);
00397 else
00398 {
00399 A.mult(x, workVector);
00400 x.update(workVector);
00401 opLU.solve(x, &y);
00402 }
00403 break;
00404 case 2:
00405 B.mult(x, y);
00406 break;
00407 default:
00408 break;
00409 }
00410 }
00411 }
00412
00413 long rvec = 1;
00414
00415
00416 const int howMnyLength = 1;
00417 char howMny = 'A';
00418
00419
00420
00421 vnl_vector<long> select(numberLanczosVecsL);
00422
00423
00424 nvalues = nconv;
00425 values = new double[nvalues];
00426 vectors = new vnl_vector<double>[nvalues];
00427
00428
00429 double *Z = new double[nvalues * matSize];
00430
00431 v3p_netlib_dseupd_(&rvec, &howMny, select.data_block(), values, Z, &matSize, &sigma, &genEigProblem,
00432 &matSize, which, &nEVL, &tolerance, resid, &numberLanczosVecsL, &V[1], &matSize, &iParam[1],
00433 &iPntr[1], &workd[1], &workl[1], &lworkl, &info,
00434 howMnyLength, genEigProblemLength, whichLength);
00435
00436
00437 int evIx;
00438 for (evIx = 0; evIx < nvalues; evIx++)
00439 {
00440 vnl_vector_ref<double> tempEVec(matSize, &Z[evIx * matSize]);
00441 vectors[evIx] = tempEVec;
00442 }
00443
00444
00445 delete[] Z;
00446 delete[] resid;
00447 delete[] V;
00448 delete[] workd;
00449 delete[] workl;
00450
00451 return info;
00452 }
00453
00454
00455
00456
00457 int vnl_sparse_symmetric_eigensystem::CalculateProduct(int n, int m,
00458 const double* p,
00459 double* q)
00460 {
00461
00462 mat->mult(n,m,p,q);
00463
00464 return 0;
00465 }
00466
00467
00468
00469 int vnl_sparse_symmetric_eigensystem::SaveVectors(int n, int m,
00470 const double* q,
00471 int base)
00472 {
00473
00474
00475 if (base == 0) {
00476 for (unsigned i=0; i<temp_store.size(); ++i)
00477 delete temp_store[i];
00478 temp_store.clear();
00479 }
00480
00481 double* temp = new double[n*m];
00482 vcl_memcpy(temp,q,n*m*sizeof(double));
00483 #ifdef DEBUG
00484 vcl_cout << "Save vectors " << base << ' ' << temp << '\n';
00485 #endif
00486
00487 temp_store.push_back(temp);
00488
00489 return 0;
00490 }
00491
00492
00493
00494 int vnl_sparse_symmetric_eigensystem::RestoreVectors(int n, int m,
00495 double* q,
00496 int base)
00497 {
00498
00499
00500 static int read_idx = 0;
00501 if (base == 0)
00502 read_idx = 0;
00503
00504 double* temp = temp_store[read_idx];
00505 vcl_memcpy(q,temp,n*m*sizeof(double));
00506 #ifdef DEBUG
00507 vcl_cout << "Restore vectors " << base << ' ' << temp << '\n';
00508 #endif
00509
00510 read_idx++;
00511 return 0;
00512 }
00513
00514
00515
00516 vnl_vector<double> vnl_sparse_symmetric_eigensystem::get_eigenvector(int i) const
00517 {
00518 assert(i>=0 && i<nvalues);
00519 return vectors[i];
00520 }
00521
00522 double vnl_sparse_symmetric_eigensystem::get_eigenvalue(int i) const
00523 {
00524 assert(i>=0 && i<nvalues);
00525 return values[i];
00526 }