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