1 /*
2 ARPACK++ v1.2 2/20/2000
3 c++ interface to ARPACK code.
4
5 MODULE ARUSPen.h.
6 Arpack++ class ARumSymPencil definition.
7
8 Modified to work with Umfpack v5.??
9 Martin Reuter
10 Date 02/28/2013
11
12 Arpack++ Author:
13 Francisco Gomes
14
15
16 ARPACK Authors
17 Richard Lehoucq
18 Danny Sorensen
19 Chao Yang
20 Dept. of Computational & Applied Mathematics
21 Rice University
22 Houston, Texas
23 */
24
25 #ifndef ARUSPEN_H
26 #define ARUSPEN_H
27
28 //#include "arch.h"
29 //#include "arerror.h"
30 //#include "lapackc.h"
31 #include "arusmat.h"
32 #include "blas1c.h"
33
34
35 template<class ARTYPE>
36 class ARumSymPencil
37 {
38
39 protected:
40
41 ARumSymMatrix<ARTYPE>* A;
42 ARumSymMatrix<ARTYPE>* B;
43 //ARumSymMatrix<ARTYPE> AsB;
44 void* Numeric;
45 int* Ap;
46 int* Ai;
47 ARTYPE* Ax;
48
49 virtual void Copy(const ARumSymPencil& other);
50
51 // void SparseSaxpy(ARTYPE a, ARTYPE x[], int xind[], int nx, ARTYPE y[],
52 // int yind[], int ny, ARTYPE z[], int zind[], int& nz);
53
54 void ExpandAsB(ARTYPE sigma);
55
56 // void SubtractAsB(ARTYPE sigma);
57 void ClearMem();
58
59 public:
60
IsFactored()61 bool IsFactored() { return (Numeric != NULL); }
62
63 void FactorAsB(ARTYPE sigma);
64
MultAv(ARTYPE * v,ARTYPE * w)65 void MultAv(ARTYPE* v, ARTYPE* w) { A->MultMv(v,w); }
66
MultBv(ARTYPE * v,ARTYPE * w)67 void MultBv(ARTYPE* v, ARTYPE* w) { B->MultMv(v,w); }
68
69 void MultInvBAv(ARTYPE* v, ARTYPE* w);
70
71 //void MultInvAsBv(ARTYPE* v, ARTYPE* w) { AsB.MultInvv(v,w); }
72 void MultInvAsBv(ARTYPE* v, ARTYPE* w);
73
74 void DefineMatrices(ARumSymMatrix<ARTYPE>& Ap, ARumSymMatrix<ARTYPE>& Bp);
75
76 //ARumSymPencil() { AsB.factored = false; }
ARumSymPencil()77 ARumSymPencil() { Numeric = NULL; Ap = NULL; Ai = NULL; Ax = NULL; }
78 // Short constructor that does nothing.
79
80 ARumSymPencil(ARumSymMatrix<ARTYPE>& Ap, ARumSymMatrix<ARTYPE>& Bp);
81 // Long constructor.
82
ARumSymPencil(const ARumSymPencil & other)83 ARumSymPencil(const ARumSymPencil& other) { Copy(other); }
84 // Copy constructor.
85
~ARumSymPencil()86 virtual ~ARumSymPencil() { }
87 // Destructor.
88
89 ARumSymPencil& operator=(const ARumSymPencil& other);
90 // Assignment operator.
91
92 };
93
94 // ------------------------------------------------------------------------ //
95 // ARumSymPencil member functions definition. //
96 // ------------------------------------------------------------------------ //
97
98
99 template<class ARTYPE>
ClearMem()100 inline void ARumSymPencil<ARTYPE>::ClearMem()
101 {
102
103 if (Numeric) umfpack_di_free_numeric (&Numeric);
104 if (Ai) delete [] Ai;
105 Ai = NULL;
106 if (Ap) delete [] Ap;
107 Ap = NULL;
108 if (Ax) delete [] Ax;
109 Ax = NULL;
110
111 } // ClearMem.
112
113
114
115 template<class ARTYPE>
Copy(const ARumSymPencil<ARTYPE> & other)116 inline void ARumSymPencil<ARTYPE>::Copy(const ARumSymPencil<ARTYPE>& other)
117 {
118 ClearMem();
119 A = other.A;
120 B = other.B;
121 // AsB = other.AsB;
122
123 } // Copy.
124
125
126 /*template<class ARTYPE>
127 void ARumSymPencil<ARTYPE>::
128 SparseSaxpy(ARTYPE a, ARTYPE x[], int xind[], int nx, ARTYPE y[],
129 int yind[], int ny, ARTYPE z[], int zind[], int& nz)
130 // A strongly sequential (and inefficient) sparse saxpy algorithm.
131 {
132
133 int ix, iy;
134
135 nz = 0;
136 if ((nx == 0) || (a == (ARTYPE)0)) {
137 copy(ny,y,1,z,1);
138 for (iy=0; iy!=ny; iy++) zind[iy] = yind[iy];
139 nz = ny;
140 return;
141 }
142 if (ny == 0) {
143 copy(nx,x,1,z,1);
144 scal(nx,a,z,1);
145 for (ix=0; ix!=nx; ix++) zind[ix] = xind[ix];
146 nz = nx;
147 return;
148 }
149 ix = 0;
150 iy = 0;
151 while (true) {
152 if (xind[ix] == yind[iy]) {
153 zind[nz] = xind[ix];
154 z[nz++] = a*x[ix++]+y[iy++];
155 if ((ix == nx)||(iy == ny)) break;
156 }
157 else if (xind[ix] < yind[iy]) {
158 zind[nz] = xind[ix];
159 z[nz++] = a*x[ix++];
160 if (ix == nx) break;
161 }
162 else {
163 zind[nz] = yind[iy];
164 z[nz++] = y[iy++];
165 if (iy == ny) break;
166 }
167 }
168 while (iy < ny) {
169 zind[nz] = yind[iy];
170 z[nz++] = y[iy++];
171 }
172 while (ix < nx) {
173 zind[nz] = xind[ix];
174 z[nz++] = x[ix++];
175 }
176
177 } // SparseSaxpy.
178
179
180 template<class ARTYPE>
181 void ARumSymPencil<ARTYPE>::ExpandAsB()
182 {
183
184 int i, j, k, n;
185 int *pcol, *irow, *index, *pos;
186 ARTYPE *value;
187
188 // Initializing variables.
189
190 n = AsB.n;
191 index = AsB.index;
192 value = AsB.value;
193 irow = &index[n+1];
194 pcol = new int[AsB.n+1];
195 pos = new int[AsB.n+1];
196 for (i=0; i<=n; i++) pcol[i] = index[i];
197 for (i=0; i<=n; i++) pos[i] = 0;
198
199 // Counting the elements in each column of AsB.
200
201 if (AsB.uplo == 'U') {
202
203 for (i=0; i!=n; i++) {
204 k = pcol[i+1];
205 if ((k!=pcol[i])&&(irow[k-1]==i)) k--;
206 for (j=pcol[i]; j<k; j++) pos[irow[j]]++;
207 }
208
209 }
210 else { // uplo == 'L'
211
212 for (i=0; i!=n; i++) {
213 k = pcol[i];
214 if ((k!=pcol[i+1])&&(irow[k]==i)) k++;
215 for (j=k; j<pcol[i+1]; j++) pos[irow[j]]++;
216 }
217
218 }
219
220 // Summing up index elements.
221
222 for (i=0; i<n; i++) pos[i+1] += pos[i];
223 for (i=n; i>0; i--) index[i] += pos[i-1];
224
225 // Expanding A.
226
227 if (AsB.uplo == 'U') {
228
229 for (i=n-1; i>=0; i--) {
230 pos[i] = index[i]+pcol[i+1]-pcol[i];
231 k = pos[i]-1;
232 for (j=pcol[i+1]-1; j>=pcol[i]; j--) {
233 value[k] = value[j];
234 irow[k--] = irow[j];
235 }
236 }
237 for (i=1; i<n; i++) {
238 k = index[i]+pcol[i+1]-pcol[i];
239 if ((k>index[i])&&(irow[k-1]==i)) k--;
240 for (j=index[i]; j<k; j++) {
241 value[pos[irow[j]]] = value[j];
242 irow[pos[irow[j]]++] = i;
243 }
244 }
245
246 }
247 else { // uplo == 'L'
248
249 for (i=n-1; i>=0; i--) {
250 k = index[i+1]-1;
251 for (j=pcol[i+1]-1; j>=pcol[i]; j--) {
252 value[k] = value[j];
253 irow[k--] = irow[j];
254 }
255 pos[i] = index[i];
256 }
257 for (i=0; i<(n-1); i++) {
258 k = index[i+1]-pcol[i+1]+pcol[i];
259 if ((k<index[i+1])&&(irow[k]==i)) k++;
260 for (j=k; j<index[i+1]; j++) {
261 value[pos[irow[j]]] = value[j];
262 irow[pos[irow[j]]++] = i;
263 }
264 }
265
266 }
267
268 AsB.nnz = index[n];
269
270 // Deleting temporary vectors.
271
272 delete[] pcol;
273 delete[] pos;
274
275 } // ExpandAsB.
276
277
278 template<class ARTYPE>
279 void ARumSymPencil<ARTYPE>::SubtractAsB(ARTYPE sigma)
280 {
281
282 int i, acol, bcol, asbcol, scol;
283
284 // Quitting function if A->uplo is not equal to B->uplo.
285
286 if ((A->uplo != B->uplo)&&(sigma != (ARTYPE)0)) {
287 throw ArpackError(ArpackError::DIFFERENT_TRIANGLES,
288 "ARumSymPencil::SubtractAsB");
289 }
290 AsB.uplo = A->uplo;
291
292 // Subtracting sigma*B from A.
293
294 AsB.index[0] = 0;
295 asbcol = 0;
296
297 for (i=0; i!=AsB.n; i++) {
298 bcol = B->pcol[i];
299 acol = A->pcol[i];
300 SparseSaxpy(-sigma, &B->a[bcol], &B->irow[bcol], B->pcol[i+1]-bcol,
301 &A->a[acol], &A->irow[acol], A->pcol[i+1]-acol,
302 &AsB.value[asbcol], &AsB.index[asbcol+AsB.n+1], scol);
303 asbcol += scol;
304 AsB.index[i+1] = asbcol;
305 }
306
307 // Expanding AsB.
308
309 ExpandAsB();
310
311 // Adding one to all elements of vector index
312 // because the decomposition function was written in FORTRAN.
313
314 for (i=0; i<=AsB.n+AsB.nnz; i++) AsB.index[i]++;
315
316 } // SubtractAsB. */
317
318
319 template<class ARTYPE>
ExpandAsB(ARTYPE sigma)320 void ARumSymPencil<ARTYPE>::ExpandAsB(ARTYPE sigma)
321 {
322 std::cout <<"ARumSymPencil::ExpandAsB(" << sigma << ") ..." << std::flush;
323
324 ClearMem();
325
326 int mynnz = 2*A->nnz+2*B->nnz;
327 if (sigma == 0.0)
328 mynnz = 2*A->nnz;
329
330 // create triples (i,j,value)
331 int * tripi = new int[mynnz];
332 int * tripj = new int[mynnz];
333 ARTYPE* tripx = new ARTYPE[mynnz];
334 if (tripi == NULL || tripj == NULL || tripx ==NULL)
335 throw ArpackError(ArpackError::PARAMETER_ERROR, "ARumSymPencil::ExpandAsB out of memory (1)");
336
337 int count = 0;
338 int i,j;
339 for (i=0; i < A->n; i++)
340 {
341 // create triplets from A
342 for (j=A->pcol[i]; j<(A->pcol[i+1]); j++)
343 {
344 tripi[count] = i;
345 tripj[count] = A->irow[j];
346 tripx[count] = A->a[j];
347 count++;
348 if (i != A->irow[j]) // not on diag
349 {
350 tripj[count] = i;
351 tripi[count] = A->irow[j];
352 tripx[count] = A->a[j];
353 count++;
354 }
355 }
356
357 if (sigma != 0.0)
358 {
359 // create triplets from -sigma B
360 for (j=B->pcol[i]; j<(B->pcol[i+1]); j++)
361 {
362 tripi[count] = i;
363 tripj[count] = B->irow[j];
364 tripx[count] = -sigma * B->a[j];
365 count++;
366 if (i != B->irow[j]) // not on diag
367 {
368 tripj[count] = i;
369 tripi[count] = B->irow[j];
370 tripx[count] = tripx[count-1];
371 count++;
372 }
373 }
374 }
375
376 }
377
378 //Write_Triplet_Matrix("A-aruspen.asc",tripi,tripj,tripx,count);
379
380 std::cout<< " ( N = " << A->n << " NNZ = " << count << " )" << std::flush;
381 //std::cout<< " size double " << sizeof(double) << " size ARTYPE " << sizeof(ARTYPE) << std::endl;
382 // convert triples (A-sigma B) to Ax Ap Ai
383 Ap = new int[A->n + 1];
384 Ai = new int[count];
385 Ax = new ARTYPE[count];
386 if (!Ap || !Ai || !Ax )
387 throw ArpackError(ArpackError::PARAMETER_ERROR, "ARumSymPencil::ExpandAsB out of memory (2)");
388
389 int status = umfpack_di_triplet_to_col (A->n, A->n, count, tripi, tripj, tripx, Ap, Ai, Ax, (int *)NULL) ;
390 if (status != UMFPACK_OK)
391 throw ArpackError(ArpackError::PARAMETER_ERROR, "ARumSymPencil::ExpandAsB triplet to col");
392
393 // cleanup
394 delete [] tripi;
395 delete [] tripj;
396 delete [] tripx;
397
398 //std::cout << std::endl << std::endl;
399 //double Control [UMFPACK_CONTROL];
400 //Control [UMFPACK_PRL] = 3;
401 //status = umfpack_di_report_matrix(A->n, A->n,Ap, Ai, Ax,0,Control);
402 //std::cout << " status: " << status << std::endl;
403 //std::cout << std::endl << std::endl;
404
405 std::cout <<" done!" << std::endl;
406 }
407
408 template<class ARTYPE>
FactorAsB(ARTYPE sigma)409 void ARumSymPencil<ARTYPE>::FactorAsB(ARTYPE sigma)
410 {
411
412 // Quitting the function if A and B were not defined.
413
414 if (!(A->IsDefined()&&B->IsDefined())) {
415 throw ArpackError(ArpackError::DATA_UNDEFINED,
416 "ARumSymPencil::FactorAsB");
417 }
418
419
420 // Subtracting sigma*B from A and storing the result
421 ExpandAsB(sigma);
422
423 // Decomposing AsB.
424 double Info [UMFPACK_INFO], Control [UMFPACK_CONTROL];
425 umfpack_di_defaults (Control) ;
426 //std::cout <<" loaded defaults" << std::endl;
427 void *Symbolic ;
428 int status = umfpack_di_symbolic (A->n, A->n, Ap, Ai, Ax, &Symbolic, Control, Info) ;
429 std::cout << " symbolic status: " << status << std::endl;
430 if (status != UMFPACK_OK)
431 throw ArpackError(ArpackError::PARAMETER_ERROR, "ARumSymPencil::FactorAsB symbolic");
432 status = umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric, Control, Info) ;
433 std::cout << " numeric status: " << status << std::endl;
434 if (status == 1)
435 {
436 std::cout << " WARNING: MATRIX IS SINGULAR " << std::endl;
437 //throw ArpackError(ArpackError::PARAMETER_ERROR, "ARumSymPencil::FactorAsB numeric (matrix singular)");
438 }
439 if (status < UMFPACK_OK)
440 {
441 std::cout << " ERROR CODE: " << status << std::endl;
442 throw ArpackError(ArpackError::PARAMETER_ERROR, "ARumSymPencil::FactorAsB numeric");
443 }
444 umfpack_di_free_symbolic (&Symbolic) ;
445
446 //exit(0);
447
448 // Decomposing AsB.
449
450 //um2fa(AsB.n, AsB.index[AsB.n], 0, false, AsB.lvalue, AsB.lindex, AsB.value,
451 // AsB.index, AsB.keep, AsB.cntl, AsB.icntl, AsB.info, AsB.rinfo);
452
453 // Handling errors.
454
455 // AsB.ThrowError();
456
457 // AsB.factored = true;
458
459 } // FactorAsB (ARTYPE shift).
460
461
462 template<class ARTYPE>
MultInvBAv(ARTYPE * v,ARTYPE * w)463 void ARumSymPencil<ARTYPE>::MultInvBAv(ARTYPE* v, ARTYPE* w)
464 {
465
466 if (!B->IsFactored()) B->FactorA();
467
468 A->MultMv(v, w);
469 copy(A->ncols(), w, 1, v, 1);
470 B->MultInvv(w, w);
471
472 } // MultInvBAv.
473
474 template<class ARTYPE>
MultInvAsBv(ARTYPE * v,ARTYPE * w)475 void ARumSymPencil<ARTYPE>::MultInvAsBv(ARTYPE* v, ARTYPE* w)
476 {
477 if (!Numeric) {
478 throw ArpackError(ArpackError::NOT_FACTORED_MATRIX,
479 "ARchSymPencil::MultInvAsBv");
480 }
481
482 // Solving A.w = v (or AsI.w = v).
483 int status = umfpack_di_solve (UMFPACK_A, Ap, Ai, Ax, w, v, Numeric, NULL, NULL) ;
484 if (status == 1)
485 {
486 std::cout << " WARNING: MATRIX IS SINGULAR " << std::endl;
487 //throw ArpackError(ArpackError::PARAMETER_ERROR, "ARumSymPencil::FactorAsB numeric (matrix singular)");
488 }
489 if (status < UMFPACK_OK)
490 {
491 std::cout << " ERROR CODE: " << status << std::endl;
492 throw ArpackError(ArpackError::PARAMETER_ERROR, "ARumSymPencil::MultInvAsBv");
493
494 }
495
496 } // MultInvAsBv
497
498 template<class ARTYPE>
499 inline void ARumSymPencil<ARTYPE>::
DefineMatrices(ARumSymMatrix<ARTYPE> & Ap,ARumSymMatrix<ARTYPE> & Bp)500 DefineMatrices(ARumSymMatrix<ARTYPE>& Ap, ARumSymMatrix<ARTYPE>& Bp)
501 {
502
503 A = &Ap;
504 B = &Bp;
505
506 if (A->n != B->n) {
507 throw ArpackError(ArpackError::INCOMPATIBLE_SIZES,
508 "ARumSymMatrix::DefineMatrices");
509 }
510
511 } // DefineMatrices.
512
513
514 template<class ARTYPE>
515 inline ARumSymPencil<ARTYPE>::
ARumSymPencil(ARumSymMatrix<ARTYPE> & Ap,ARumSymMatrix<ARTYPE> & Bp)516 ARumSymPencil(ARumSymMatrix<ARTYPE>& Ap, ARumSymMatrix<ARTYPE>& Bp)
517 {
518 Numeric = NULL;
519 Ap = NULL;
520 Ai = NULL;
521 Ax = NULL;
522
523 //AsB.factored = false;
524 DefineMatrices(Ap, Bp);
525
526
527 } // Long constructor.
528
529
530 template<class ARTYPE>
531 ARumSymPencil<ARTYPE>& ARumSymPencil<ARTYPE>::
532 operator=(const ARumSymPencil<ARTYPE>& other)
533 {
534
535 if (this != &other) { // Stroustrup suggestion.
536 Copy(other);
537 }
538 return *this;
539
540 } // operator=.
541
542
543 #endif // ARUSPEN_H
544