source: sundry/stable/1.1/sundry/SMVkernel.h @ 5999

Revision 5999, 35.0 KB checked in by wehart, 5 years ago (diff)

Merged revisions 5281-5998 via svnmerge from
 https://software.sandia.gov/svn/public/acro/sundry/trunk

........

r5296 | jberry | 2007-12-19 10:22:19 -0700 (Wed, 19 Dec 2007) | 2 lines


Use the multipleRandomizedSensorPlacements correctly.

........

r5376 | lafisk | 2008-05-23 13:56:37 -0600 (Fri, 23 May 2008) | 2 lines


Add link flags for static executables.

........

r5834 | wehart | 2008-11-23 20:01:45 -0700 (Sun, 23 Nov 2008) | 2 lines


Update of Copyright assertions.

........

r5853 | wehart | 2008-11-24 08:42:59 -0700 (Mon, 24 Nov 2008) | 2 lines


Update to change the license on these files to CPL.

........

r5901 | wehart | 2008-12-01 20:32:01 -0700 (Mon, 01 Dec 2008) | 2 lines


Portability fixes for gcc 4.3.1

........

  • Property svn:executable set to *
Line 
1/*  _________________________________________________________________________
2 *
3 *  Acro: A Common Repository for Optimizers
4 *  Copyright (c) 2008 Sandia Corporation.
5 *  This software is distributed under the BSD License.
6 *  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
7 *  the U.S. Government retains certain rights in this software.
8 *  For more information, see the README.txt file in the top Acro directory.
9 *  _________________________________________________________________________
10 */
11
12//
13//  SMVkernel.h
14//        Description: Thread-safe sparse matrix data structures
15//                     and algorithms
16//        Author:  Robert Heaphy
17//        Date:    2005
18//
19#ifndef __SMVKERNEL_H__
20#define __SMVKERNEL_H__
21
22#include <stdio.h>
23#include <stdlib.h>
24#include <math.h>
25#include <malloc.h>
26
27//#include <mtgl/graph.h>
28
29#ifdef __MTA__
30 #define BRAKET
31 #include <sys/mta_task.h>
32 #include <machine/runtime.h>
33#else
34 #define BRAKET <T>
35#endif
36
37// Current ISSUE: error handling is erratic
38// Current ISSUE: look for NOT IMPLEMENTED and implement
39
40
41// Base classes-forward declarations
42template <class T> class VectorBase;
43template <class T> class MatrixBase;
44
45// Derived classes-forward declarations
46template <class T> class DenseVector;
47
48template <class T> class SparseMatrixCSR;
49template <class T> class SparseMatrixCSC;
50template <class T> class SparseMatrixCOO;
51
52template <class T> DenseVector<T>
53operator* (const SparseMatrixCSR<T>& a, const DenseVector<T>& b);
54template <class T> DenseVector<T> diagonal (const SparseMatrixCSR<T>& a);
55template <class T> DenseVector<T> Transpose_SMVm (const SparseMatrixCSR<T>&,
56                                                  const DenseVector<T>&);
57template <class T>
58DenseVector<T>
59operator* (const SparseMatrixCSC<T>&, const DenseVector<T>&);
60     
61template <class T>
62DenseVector<T>
63operator* (const SparseMatrixCOO<T>&, const DenseVector<T>&);
64     
65template <class T>
66DenseVector<T>
67operator* (int const, const DenseVector<T>&);
68   
69template <class T>
70DenseVector<T> operator* (double const, const DenseVector<T>&);
71
72template <class T>
73SparseMatrixCSR<T> operator*(int const, const SparseMatrixCSR<T>&);   
74
75template <class T>
76SparseMatrixCSR<T>
77     operator* (double const, const SparseMatrixCSR<T>&);
78   
79template <class T>
80T
81operator* (const DenseVector<T>&, const DenseVector<T>&);
82     
83template <class T>
84DenseVector<T>
85diagonal (const SparseMatrixCSR<T>&);
86   
87template <class T>
88DenseVector<T>
89Transpose_SMVm (const SparseMatrixCSR<T>&, const DenseVector<T>&);
90
91/***********************  Base Classes  ****************************/
92
93template <class T> class VectorBase
94{
95  protected:
96    int length;               
97    T* values;                                         // vector elements
98   
99  public: 
100    VectorBase() : length(0), values(0)  {
101    }
102     
103    VectorBase (int const n) : length(n)  {
104      this->values = (T*) malloc (n * sizeof(T));            // values(new T[n])
105     
106    T*  const this_values = this->values;         // force MTA to place on stack
107    T   const zero        = T();                  // force MTA to place on stack
108    int const finish      = this->length;         // force MTA to place on stack
109    #pragma mta assert parallel
110    for (int i = 0; i < finish; i++)
111      this_values[i] = zero;       
112    }
113   
114    VectorBase (const VectorBase<T>& a) : length(a.length)  {
115      this->values = (T*) malloc (a.length * sizeof(T)); //values(new T[a.length])
116     
117      int const stop        = a.length;            // force MTA to place on stack
118      T*  const this_values = this->values;        // force MTA to place on stack
119      T*  const a_values    = a.values;            // force MTA to place on stack
120      #pragma mta assert parallel 
121      for (int i = 0; i < stop; i++)
122        this_values[i] = a_values[i];
123    }
124     
125    ~VectorBase()  {
126      if (this->values) free (this->values);                  // delete[] values
127      this->values = 0;
128    }
129
130    void clear()
131    {
132        if(this->values) free(this->values); this->values = NULL;
133    }
134         
135    void VectorPrint (char const* name) const;
136};
137
138
139
140template <class T> class MatrixBase
141{
142  protected:
143    int nRow, nCol, nNonZero;     // number of rows, columns, and non zeros
144    T*   values;                  // non zero matrix elements
145    int* index;                   // generally used for sparse storage index
146    bool data_owned;
147   
148  public:
149    MatrixBase() : nRow(0), nCol(0), nNonZero(0), values(0), index(0),
150                   data_owned(true)  { }
151     
152    MatrixBase(int const row, int const col, int const count, bool own=true)
153     : nRow(row), nCol(col), nNonZero(count), index(0), data_owned(own) 
154    {
155        //printf("MatrixBase()\n");
156        values = (T*) malloc (count * sizeof(T));  // values(new T[count])
157        //printf("MatrixBase() alloc'd\n");
158
159        T*  const this_values = this->values;      // force MTA to place on stack
160        T   const zero        = T();
161        int const finish      = this->nNonZero;
162       
163        #pragma mta assert parallel
164        for (int i = 0; i < finish; i++)
165        {
166            this_values[i] = i;
167        }
168    }
169     
170    ~MatrixBase() {
171        if(data_owned && this->values) free(this->values);
172        this->values = 0;
173        if(data_owned && this->index)  free(this->index); 
174        this->index  = 0;
175    }
176
177
178    MatrixBase(const MatrixBase<T>& a)
179        : nRow(a.nRow),
180          nCol(a.nCol),
181          nNonZero(a.nNonZero),
182          /* values(new T[a.nNonZero]),*/ 
183          index(0),
184          data_owned(true)  // not debugged!
185    {
186        this->values = (T*) malloc (a.nNonZero * sizeof (T));
187        T*  const this_values = this->values;  // force MTA to place on stack
188        T*  const a_values    = a.values;      // force MTA to place on stack
189        int const stop        = a.nNonZero;
190       
191        #pragma mta assert parallel 
192        for (int i = 0; i < stop; i++)
193        {
194            this_values[i] = a_values[i];
195        }
196    }
197                 
198    void MatrixPrint (char const* name);
199};
200
201
202
203/***************************  DenseVector *************************************/
204
205template <class T = double> class DenseVector : public VectorBase<T>
206{
207  public:   
208    DenseVector()  {}
209     
210    DenseVector (int const n) : VectorBase<T> (n)  {}
211
212    DenseVector (const DenseVector<T>& a) : VectorBase<T> (a)  {}
213     
214    ~DenseVector()  {}
215
216    DenseVector<T>& operator= (const DenseVector<T>& a) 
217    {
218        if (this->length != a.length)
219            printf("DenseVector copy assignment error\n");
220       
221        if (this != &a)
222        {
223            int const stop        = this->length;  // force MTA to place on stack
224            T*  const this_values = this->values;  // force MTA to place on stack
225            T*  const a_values    = a.values;      // force MTA to place on stack
226       
227            #pragma mta assert parallel 
228            for (int i = 0; i < stop; i++)
229            {
230                this_values[i] = a_values[i];
231            }
232        }
233        return *this;
234    }
235
236
237    const T operator[](int idx)
238    {
239        if(idx<0 || idx>this->length)
240        {
241            fprintf(stderr, "error: index %d out of bounds in DenseVector\n",idx);
242            exit(1);
243        }
244        return( this->values[idx] );
245    }
246   
247     
248    DenseVector<T>& operator+= (const DenseVector<T>& a) 
249    {
250        int const stop        = this->length;    // force MTA to place on stack
251        T*  const this_values = this->values;    // force MTA to place on stack
252        T*  const a_values    = a.values;        // force MTA to place on stack
253        #pragma mta assert parallel 
254        for (int i = 0; i < stop; i++)
255        {
256            // safe_incr(this_values[i], a_values[i]);
257            this_values[i] += a_values[i];
258        }
259        return *this;
260    }
261
262    DenseVector<T>& operator-= (const DenseVector<T>& a) 
263    {
264        int const stop        = this->length;    // force MTA to place on stack
265        T*  const this_values = this->values;    // force MTA to place on stack
266        T*  const a_values    = a.values;        // force MTA to place on stack
267       
268        #pragma mta assert parallel 
269        for (int i = 0; i < stop; i++)
270        {
271            // safe_incr(this->values[i], -a.values[i]);
272            this->values[i] += -a.values[i];
273        }
274        return *this;
275    }
276                   
277    DenseVector<T>& operator() () 
278    {
279        int const stop        = this->length;    // force MTA to place on stack
280        T*  const this_values = this->values;    // force MTA to place on stack
281        T   const zero        = T();
282       
283        #pragma mta assert parallel 
284        for (int i = 0; i < stop; i++)
285        {
286            this_values[i] = zero;
287        }
288        return *this;
289    }
290
291    T norm2 () const 
292    {
293        T   temp              = T();
294        int const stop        = this->length;
295        T*  const this_values = this->values;    // force MTA to place on stack
296       
297        #pragma mta assert parallel
298        for (int i = 0 ; i < stop; i++)
299        {
300            // safe_incr(temp, this_values[i] * this_values[i]);
301            temp += this_values[i]*this_values[i];
302        }
303        return (T) sqrt(temp);
304    }
305     
306    T norm_inf () const
307    {
308        T temp                = T();
309        int const stop        = this->length;
310        T*  const this_values = this->values;   // force MTA to place on stack
311       
312        #pragma mta assert parallel
313        for (int i = 0; i < stop; i++)
314        {
315            if (fabs(this_values[i]) > temp)
316            {
317                temp = (T) fabs(this_values[i]);
318            }
319        }
320        return temp;
321    }
322     
323    int VectorLength() const 
324    {
325        return this->length;
326    } 
327     
328    friend DenseVector<T>
329     operator* BRAKET(const SparseMatrixCSR<T>&, const DenseVector<T>&);
330     
331    DenseVector<T> friend
332     operator* BRAKET (const SparseMatrixCSC<T>&, const DenseVector<T>&);
333     
334    DenseVector<T> friend
335     operator* BRAKET (const SparseMatrixCOO<T>&, const DenseVector<T>&);
336     
337    DenseVector<T> friend
338     operator* BRAKET (int const, const DenseVector<T>&);
339   
340    DenseVector<T> friend
341     operator* BRAKET (double const, const DenseVector<T>&);
342   
343    T friend
344     operator* BRAKET (const DenseVector<T>&, const DenseVector<T>&);
345     
346    DenseVector<T> friend
347     diagonal BRAKET (const SparseMatrixCSR<T>&);
348   
349    DenseVector<T> friend
350     Transpose_SMVm BRAKET (const SparseMatrixCSR<T>&, const DenseVector<T>&);
351   
352    // approximate solver, asolve, derived from  asolve.c found with linbcg.c
353    // in "Numerical Recipes in C", second edition, Press, Vetterling, Teukolsky,
354    // Flannery, pp 86-89
355    DenseVector<T> asolve (DenseVector<T>& a)  {
356      DenseVector<T> temp(this->length);
357 
358      int const stop        = this->length;
359      T*  const temp_values = temp.values;   // force MTA to place on stack
360      T*  const this_values = this->values;  // force MTA to place on stack
361      T*  const a_values    = a.values;      // force MTA to place on stack
362      T   const zero        = T();
363      #pragma mta assert parallel
364      for (int i = 0; i < stop; i++)
365        temp_values[i] = (this_values[i] != zero)
366         ? a_values[i]/this_values[i] : a_values[i];
367      return temp;
368    }
369         
370    // the following methods are for development and may be removed/modified
371    /* DenseVector::fill() */
372    void fill (T const* const val)  {
373      int const stop        = this->length;
374      T*  const this_values = this->values;  // to force MTA to place on stack
375      #pragma mta assert parallel
376      for (int i = 0; i < stop; i++) {
377        this_values[i] = val[i];
378      }
379    }
380     
381    void VectorPrint(char const* name) 
382    {
383        printf("DenseVector Print %s length %d\n", name,this->length);
384     
385        int minimum = (this->length < 10) ? this->length : 10;
386     
387        printf("DenseVector Values: ");
388        for (int i = 0; i < minimum; i++)
389        {
390            printf("%g, ", this->values[i]);
391        }
392        printf("\n");
393    }
394
395};
396
397
398
399template <class T> DenseVector<T> inline
400operator* (int const a, const DenseVector<T>& b) 
401{
402    DenseVector<T> temp (b.length);
403
404    int const stop        = b.length;
405    T*  const temp_values = temp.values;       // force MTA to place on stack
406    T*  const b_values    = b.values;          // force MTA to place on stack
407
408    #pragma mta assert parallel
409    for (int i = 0; i < stop; i++)
410    {
411        temp_values[i] = a * b_values[i];
412    }
413    return temp;
414}
415
416template <class T> DenseVector<T> inline
417operator* (double const a, const DenseVector<T>& b) 
418{
419    DenseVector<T> temp (b.length);
420 
421    int const stop        = b.length;
422    T*  const temp_values = temp.values;      // force MTA to place on stack
423    T*  const b_values    = b.values;          // force MTA to place on stack
424   
425    #pragma mta assert parallel
426    for (int i = 0; i < stop; i++)
427    {
428        temp_values[i] = a * b_values[i];
429    }
430    return temp;
431}
432
433template <class T> T inline
434operator* (const DenseVector<T>& a, const DenseVector<T>& b) {
435   
436   T temp = T();
437   int const stop     = a.length;
438   T*  const a_values = a.values;            // force MTA to place on stack
439   T*  const b_values = b.values;            // force MTA to place on stack
440   #pragma mta assert parallel
441   for (int i = 0; i < stop; i++)
442      //safe_incr(temp, a_values[i] * b_values[i]) ;
443      temp += a_values[i]*b_values[i];
444   return temp;
445}
446     
447template <class T> DenseVector<T> 
448operator- (const DenseVector<T>& a, const DenseVector<T>& b) {
449   DenseVector<T> temp(a);
450   temp -= b;
451   return temp;
452}
453
454template <class T> DenseVector<T> 
455operator+ (const DenseVector<T>& a, const DenseVector<T>& b) {
456   DenseVector<T> temp(a);
457   temp += b;
458   return temp;
459}   
460 
461   
462   
463/***************************** CSR  SparseMatrix ******************************/
464/* CSR : Compressed Sparse Row                                                */
465template <class T = double> class SparseMatrixCSR : MatrixBase<T>
466{
467  private:
468    int* columns;
469   
470  public:
471    SparseMatrixCSR() : MatrixBase<T>(), columns(0)
472    {
473    }
474     
475    SparseMatrixCSR(int const row, int const col, int const count)
476                    : MatrixBase<T> (row, col, count,true)
477    {
478         // columns=new int[count];
479        this->columns=(int*)malloc(count *sizeof(int));
480
481        // index = new int [nRow+1];
482        this->index  =(int*)malloc((this->nRow+1)*sizeof(int));
483       
484        // force MTA to place on stack
485        int* const this_columns = this->columns;
486        int  const finish       = this->nNonZero;
487        #pragma mta assert parallel
488        for (int i = 0; i < finish; i++)
489            this_columns[i] = 0;       
490       
491        // force MTA to place on stack
492        int* const this_index = this->index;
493        int  const stop       = this->nRow+1;
494        #pragma mta assert parallel
495        for (int i = 0; i < stop; i++)
496            this_index[i] = 0;
497    }
498
499    SparseMatrixCSR(const int row, const int col, const int count,
500                    int* indx, T* val, int*  cols) 
501        : MatrixBase<T>(row, col, count,false),
502                MatrixBase<T>::index(indx), columns(cols)
503    {
504    }
505     
506    SparseMatrixCSR(const SparseMatrixCSR<T>& a) : MatrixBase<T> (a)
507    {
508        columns=(int*)malloc(this->nNonZero *sizeof (int));  //columns=new int[nNonZero];
509        this->index  =(int*)malloc((this->nRow+1) *sizeof (int));  //index  =new int[nRow+1];
510     
511        int  const stop         = this->nNonZero;
512        int* const this_columns = this->columns;   // force MTA to place on stack
513        int* const a_columns    = a.columns;       // force MTA to place on stack
514        int*   const a_index     = a.index;        // force MTA to place on stack
515       
516        #pragma mta assert parallel
517        for (int i = 0; i < stop; i++)
518        {
519            this_columns[i] = a_columns[i];
520        }
521       
522        int* const this_index = this->index;       // force MTA to place on stack
523        int  const end        = this->nRow+1;      // force MTA to place on stack
524        #pragma mta assert parallel
525        for (int i = 0; i < end; i++)
526        {
527            this_index [i] = a_index[i];
528        }
529    }
530     
531    ~SparseMatrixCSR() 
532    {
533        if (this->data_owned && this->columns) free (this->columns);
534        if (this->data_owned && this->index)   free (this->index);
535        this->columns = 0;
536        this->index   = 0;
537    }
538
539    void init(const int row, const int col, const int count,
540                    int* indx, T* vals, int*  cols) 
541    {
542        clear();
543        this->nRow = row;
544        this->nCol = col;
545        this->nNonZero = count;
546        this->index = indx;
547        this->columns = cols;
548        this->values = vals;
549        this->data_owned = false;
550    }
551
552
553    void clear() {
554        if(this->data_owned && this->values) free(this->values);
555        this->values = 0;
556        if(this->data_owned && this->index)  free(this->index); 
557        this->index  = 0;
558    }
559
560
561    SparseMatrixCSR<T>& operator= (const SparseMatrixCSR<T>& a)  {
562        if (this != &a) 
563        {
564            this->nRow     = a.nRow;
565            this->nCol     = a.nCol;
566            this->nNonZero = a.nNonZero;
567       
568//          delete[] columns;   columns = new int [this->nNonZero];
569//          delete[] index;     index   = new int [this->nRow+1];
570//          delete[] values;    values  = new T   [this->nNonZero];
571       
572            if (columns)
573                free (columns); 
574            columns = (int*) malloc(this->nNonZero * sizeof(int));
575       
576            if (this->index)   
577                free (this->index);
578            this->index = (int*) malloc((this->nRow+1) * sizeof(int));
579       
580            if (this->values) 
581                free (this->values);
582            this->values = (T*) malloc(this->nNonZero * sizeof(T));
583     
584            int  const stop         = this->nNonZero;// force MTA to place on stack
585            int* const this_columns = this->columns; // force MTA to place on stack
586            T*   const this_values  = this->values;  // force MTA to place on stack
587            int* const a_columns    = a.columns;     // force MTA to place on stack
588            T*   const a_values     = a.values;      // force MTA to place on stack
589            #pragma mta assert parallel
590            for (int i = 0; i < stop; i++) 
591            {
592                this_columns[i] = a_columns[i];
593            this_values[i]  = a_values[i];
594            }
595       
596            int  const end        = this->nRow + 1;  // force MTA to place on stack
597            int* const this_index = this->index;     // force MTA to place on stack
598            int* const a_index    = a.index;         // force MTA to place on stack
599            #pragma mta assert parallel
600            for (int i = 0; i < end; i++)
601            {
602                this_index[i] = a_index[i];
603            }
604        }
605        return *this;
606    }
607     
608    SparseMatrixCSR<T>& operator= (const SparseMatrixCSC<T>& a) 
609    {
610        printf("NOT IMPLEMENTED Converting from CSC to CSR\n");
611        return *this;
612    }
613     
614    DenseVector<T> friend
615     operator* <T> (const SparseMatrixCSR<T>&, const DenseVector<T>&);
616     
617    SparseMatrixCSR<T> friend
618     operator* <T> (double const, const SparseMatrixCSR<T>&);
619   
620    SparseMatrixCSR<T> friend
621     operator* <T> (int const, const SparseMatrixCSR<T>&);   
622     
623    DenseVector<T> friend
624     diagonal BRAKET (const SparseMatrixCSR<T>&);
625   
626    DenseVector<T> friend
627     Transpose_SMVm BRAKET (const SparseMatrixCSR<T>&, const DenseVector<T>&);
628         
629    /* SparseMatrixCSR fill()
630     *
631     */
632    void fill (int const* const indx,
633               T   const* const val,
634               int const* const cols) 
635    {
636        int  const stop       = this->nRow;     // force MTA to place on stack
637        int* const this_index = this->index;    // force MTA to place on stack
638
639        #pragma mta assert parallel
640        for (int rows=0; rows < stop; rows++)
641        {
642            this_index[rows] = indx[rows];
643        }
644       
645        this_index[this->nRow] = this->nNonZero;
646     
647        int  const end          = this->nNonZero;  // force MTA to place on stack
648        T*   const this_values  = this->values;    // force MTA to place on stack
649        int* const this_columns = this->columns;   // force MTA to place on stack
650       
651        #pragma mta assert parallel
652        for (int i=0; i < end; i++) 
653        {
654            this_values[i]  = val[i];
655            this_columns[i] = cols[i];
656        }
657    }
658   
659    T* col_values_begin(int row)
660    {
661        if ((row >= 0) && (row < this->nRow))
662                return &this->values[this->index[row]];
663        else
664                return 0;
665    }
666
667    T* col_values_end(int row)
668    {
669        int ind = row+1;
670        if ((ind >= 0) && (ind <= this->nRow))
671                return &this->values[this->index[ind]];
672        else
673                return 0;
674    }
675
676    int* col_indices_begin(int row)
677    {
678        if ((row >= 0) && (row < this->nRow))
679                return &this->columns[this->index[row]];
680        else
681                return 0;
682    }
683
684    int* col_indices_end(int row)
685    {
686        int ind = row+1;
687        if ((ind >= 0) && (ind <= this->nRow))
688                return &this->columns[this->index[ind]];
689        else
690                return 0;
691    }
692   
693    void MatrixPrint (char const* name) const 
694    {
695        printf("SparseMatrixCSR Print %s row %d col %d\n",
696                name,this->nRow,this->nCol);
697    }
698   
699   
700    int MatrixRows() const {
701      return this->nRow;
702    }
703   
704
705    int MatrixCols() const {
706      return this->nCol;
707    }     
708};
709
710template <class T> DenseVector<T>
711operator* (const SparseMatrixCSR<T>& a, const DenseVector<T>& b) {
712 
713  if (a.nCol != b.length) 
714  {
715      printf("INCOMPATIBLE SparseMatrixCSR * DenseVector multiplication\n");
716      exit(1);
717  }
718   
719  DenseVector<T> temp(b.length);
720  T*  const temp_values = temp.values;         // force MTA to place on stack
721  T   const zero        = T();
722  int const finish      = temp.length;
723  #pragma mta assert parallel
724  for (int i = 0; i < finish; i++)
725    temp_values[i] = zero;
726   
727  #ifdef __MTA__ 
728    int starttimer = mta_get_clock(0);
729  #endif
730
731  int  const stop      = a.nRow;
732  int* const a_index   = a.index;            // force MTA to place on stack
733  T*   const a_values  = a.values;           // force MTA to place on stack
734  T*   const b_values  = b.values;           // force MTA to place on stack
735  int* const a_columns = a.columns;          // force MTA to place on stack
736  #pragma mta assert parallel 
737  for (int row = 0; row < stop; row++) {
738    #pragma mta trace "next_row"
739    int const start  = a_index[row];
740    int const finish = a_index[row+1];
741 
742    for (int i = start; i < finish; i++)
743    {
744        temp_values[row] += a_values[i]*b_values[ a_columns[i] ];
745    }
746  }
747     
748#ifdef __MTA__
749int stoptimer = mta_get_clock(starttimer);
750// printf("MVm total time %g\n", stoptimer/220000000.0);
751#endif
752     
753  return temp;
754}
755
756template <class T> SparseMatrixCSR<T>
757operator* (int const a, const SparseMatrixCSR<T>& b) {     
758  SparseMatrixCSR<T> temp (b.nRow, b.nCol, b.nNonZero);
759 
760  int const stop        = b.nNonZero;
761  T*  const temp_values = temp.values;         // force MTA to place on stack
762  T*  const b_values    = b.values;            // force MTA to place on stack
763  #pragma mta assert parallel
764  for (int i = 0; i < stop; i++)
765    temp_values[i] = a * b_values[i];
766  return temp;
767}
768
769template <class T> SparseMatrixCSR<T>
770operator* (double const a, const SparseMatrixCSR<T>& b) {     
771  SparseMatrixCSR<T> temp (b.nRow, b.nCol, b.nNonZero);
772 
773  int const stop = b.nNonZero;
774  T* temp_values = temp.values;               // force MTA to place on stack
775  T* b_values    = b.values;                  // force MTA to place on stack
776  #pragma mta assert parallel
777  for (int i = 0; i < stop; i++)
778    temp_values[i] = a * b_values[i];
779  return temp;
780}
781
782template <class T> DenseVector<T> diagonal (const SparseMatrixCSR<T>& a)  {
783  if (a.nRow != a.nCol) {
784      printf("diagonal called on non square matrix\n");
785      printf("nRow=%d, nCol=%d\n", a.nRow, a.nCol);
786      exit(1);
787  }
788   
789  DenseVector<T> temp(a.nRow);
790  temp();
791 
792  int const finish = a.nRow;
793  int* const a_index   = a.index;              // force MTA to place on stack
794  int* const a_columns = a.columns;            // force MTA to place on stack
795  T* const a_values    = a.values;             // force MTA to place on stack
796  T* const temp_values = temp.values;          // force MTA to place on stack
797  #pragma mta assert parallel
798  #pragma mta loop future
799  for (int row = 0; row < finish; row++)  {
800    int const start = a_index[row];
801    int const stop  = a_index[row+1];
802   
803    #pragma mta assert parallel
804    for (int i = start; i < stop; i++)
805      if (row == a_columns[i]) {
806        temp_values[row] = a_values[i];
807        #ifndef __MTA__
808          break;
809        #endif
810      }
811  }
812  return temp;
813}
814 
815template <class T> DenseVector<T> Transpose_SMVm (const SparseMatrixCSR<T>& a,
816  const DenseVector<T>& b)  {
817     
818  if (a.nCol != b.length) 
819  {
820      printf("INCOMPATIBLE Transpose (SparseMatrixCSR) * DenseVector multiplication\n");
821      exit(1);
822  }
823   
824  DenseVector<T> temp(b.length);
825  T* const temp_values = temp.values;        // force MTA to place on stack
826  T* const a_values    = a.values;           // force MTA to place on stack
827  T* const b_values    = b.values;           // force MTA to place on stack
828  int* const a_index   = a.index;            // force MTA to place on stack
829  int* const a_columns = a.columns;          // force MTA to place on stack
830  int const stop = temp.length;
831     
832  #pragma mta assert parallel
833  T const zero = T();
834  for (int i = 0; i < stop; i++)
835    temp_values[i] = zero;
836   
837  int const finish = temp.length;
838
839  #pragma mta assert parallel   
840  for (int row = 0; row < finish; row++)  {
841    int const start = a_index[row];            // force MTA to place on stack
842    int const stop  = a_index[row+1];          // force MTA to place on stack
843   
844    for (int i = start; i < stop; i++) {
845        //mt_incr(temp_values[a_columns[i]], a_values[i] * b_values[row]);
846        temp_values[ a_columns[i] ] += a_values[i] * b_values[row];
847    }
848  }
849  return temp;         
850}
851
852
853
854/***************************** CSC  SparseMatrix ******************************/
855
856template <class T = double> class SparseMatrixCSC : MatrixBase<T> {
857  private:
858    int *rows;
859  public:
860    SparseMatrixCSC<T> () : MatrixBase<T> () {
861      rows = 0;
862    }
863       
864    SparseMatrixCSC<T> (int const row, int const col, int const count)
865     : MatrixBase<T> (row, col, count)  {
866       rows  = (int*) malloc (count    * sizeof(int)); // rows (new int [count])
867       this->index = (int*) malloc ((this->nCol+1) * sizeof(int)); // index(new int [this->nCol+1])
868    }
869       
870    SparseMatrixCSC<T> (const SparseMatrixCSC<T>& a) : MatrixBase<T> (a)  {
871      rows =(int*)malloc(this->nNonZero *sizeof(int)); // rows(new int[nNonZero])
872      this->index=(int*)malloc((this->nCol+1) *sizeof(int)); // index(new int[nCol+1])   
873     
874      int* const a_rows    = a.rows;            // force MTA to place on stack
875      int* const this_rows = this->rows;        // force MTA to place on stack
876      int  const stop      = this->nNonZero;
877      #pragma mta assert parallel
878      for (int i = 0; i < stop; i++)
879        this_rows[i] = a_rows[i];
880       
881      int* const this_index = this->index;      // force MTA to place on stack
882      int* const a_index    = a.index;          // force MTA to place on stack
883      int  const end        = this->nCol + 1;
884      #pragma mta assert parallel
885      for (int i = 0; i < end; i++)
886        this_index[i] = a_index[i];
887    }
888     
889    ~SparseMatrixCSC<T> ()  {   
890        if (this->data_owned && this->rows)  free (this->rows);
891        if (this->data_owned && this->index) free (this->index);
892        this->rows  = 0;
893        this->index = 0;
894    }
895 
896    SparseMatrixCSC<T>& operator= (const SparseMatrixCSR<T>& a) 
897    {
898        printf("NOT IMPLEMENTED: Converting from CSR to CSC\n");
899        return *this;
900    }
901     
902    SparseMatrixCSC<T>& operator= (const SparseMatrixCSC<T>& a)  {
903      if (this != &a)  {
904        this->nRow     = a.nRow;
905        this->nCol     = a.nCol;
906        this->nNonZero = a.nNonZero;
907     
908        int* const this_rows   = this->rows;       // force MTA to place on stack
909        int* const a_rows      = a.rows;           // force MTA to place on stack
910        T*   const this_values = this->values;     // force MTA to place on stack
911        T*   const a_values    = a.values;         // force MTA to place on stack
912        int  const stop        = this->nNonZero;
913        #pragma mta assert parallel
914        for (int i = 0; i < stop; i++)  {
915          this_rows[i]   = a_rows[i];
916          this_values[i] = a_values[i];
917        }
918       
919        int* const this_index = this->index;      // force MTA to place on stack
920        int* const a_index    = a.index;          // force MTA to place on stack
921        int  const end        = this->nCol + 1;
922        #pragma mta assert parallel
923        for (int i = 0; i < end; i++)
924          this_index[i] = a_index[i];
925      }
926      return *this;
927    }
928     
929    SparseMatrixCSC<T>& operator() (SparseMatrixCSR<T>&);
930   
931    friend DenseVector<T>
932     operator* <T> (const SparseMatrixCSC<T>&, const DenseVector<T>&);
933   
934    void MatrixPrint (char const* name) 
935    {
936        printf("SparseMatrixCSC Print %s row %d col %d\n",
937                name,this->nRow,this->nCol);
938    }
939     
940    int MatrixRows() const {
941      return this->nRow;
942    }
943     
944    int MatrixCols() const {
945      return this->nCol;
946    }
947};
948
949
950
951template <class T> DenseVector<T>
952operator* (const SparseMatrixCSC<T>& a, const DenseVector<T>& b) {
953  if (a.nCol != b.length) 
954  {
955      printf("INCOMPATIBLE SparseMatrixCSC * DenseVector multiplication\n");
956      exit(1);
957  }
958   
959  DenseVector<T> temp(b.length);
960  T*  const temp_values = temp.values;       // force MTA to place on stack
961  int const istop       = temp.length;
962  T   const zero        = T();
963  #pragma mta assert parallel
964  for (int i = 0; i < istop; i++)
965    temp_values[i] = zero;
966   
967  int  const colstop  = a.nCol;
968  int* const a_index  = a.index;          // force MTA to place on stack
969  int* const a_rows   = a.rows;           // force MTA to place on stack
970  T*   const a_values = a.values;         // force MTA to place on stack
971  T* const b_values    = b.values;           // force MTA to place on stack
972  #pragma mta assert parallel
973  #pragma mta loop future
974  for (int col = 0; col < colstop; col++) {
975    int const start = a_index[col];
976    int const stop  = a_index[col+1];
977 
978    #pragma mta assert parallel
979    for (int i = start; i < stop; i++) {
980        //mt_incr(temp_values[a_rows[i]], a_values[i] * b_values[a_rows[i]]);
981        temp_values[ a_rows[i] ] += a_values[i] * b_values[ a_rows[i] ];
982    }
983  }
984  return temp;
985}
986
987
988
989/***************************** COO  SparseMatrix ******************************/
990/* COO : */
991
992template <class T = double> class SparseMatrixCOO : MatrixBase<T>
993{
994  private:
995    int *columns;
996    int *rows;
997   
998  public:
999    SparseMatrixCOO() : MatrixBase<T> ()  {
1000      this->columns = 0;
1001      this->rows    = 0;
1002      }
1003     
1004    SparseMatrixCOO(int const row, int const col, int const nnz)
1005     : MatrixBase<T> (row, col, nnz)  {
1006        rows    = (int*) malloc (nnz * sizeof(int));   // rows    = new int[nnz]
1007        columns = (int*) malloc (nnz * sizeof(int));   // columns = new int[nnz]
1008    }
1009       
1010    SparseMatrixCOO (SparseMatrixCOO<T>& a) : MatrixBase<T> (a)  {
1011      rows   =(int*)malloc(a.nNonZero*sizeof(int)); //rows   =new int[a.nNonZero]
1012      columns=(int*)malloc(a.nNonZero*sizeof(int)); //columns=new int[a.nNonZero]
1013     
1014      int  const stop         = this->nNonZero;
1015      T*   const this_values  = this->values;     // MTA -to force onto stack
1016      T*   const a_values     = a.values;         // MTA -to force onto stack
1017      int* const this_rows    = this->rows;       // MTA -to force onto stack
1018      int* const this_columns = this->columns;    // MTA -to force onto stack
1019      int* const a_rows       = a.rows;           // MTA -to force onto stack
1020      int* const a_columns    = a.columns;        // MTA -to force onto stack
1021      #pragma mta assert parallel
1022      for (int i = 0; i < stop; i++)  {
1023        this_values[i]  = a_values[i];
1024        this_rows[i]    = a_rows[i];
1025        this_columns[i] = a_columns[i];
1026      }
1027    }
1028     
1029    ~SparseMatrixCOO()  {
1030       if (this->columns)  free (this->columns);         // delete[] columns
1031       if (this->index)    free (this->index);           // delete[] index
1032       if (this->rows)     free (this->rows);            // delete[] rows
1033       this->columns = 0;
1034       this->index   = 0;
1035       this->rows    = 0;
1036    }
1037     
1038    SparseMatrixCOO<T>& operator= (const SparseMatrixCOO<T>& a)  {
1039      if (this != &a)  {
1040        this->nRow     = a.nRow;
1041        this->nCol     = a.nCol;
1042        this->nNonZero = a.nNonZero;
1043        this->index    = 0;
1044       
1045        if (rows)   free(rows);                      // delete[] rows;
1046        rows=(int*)malloc(a.nNonZero*sizeof(int));   // rows = new int[a.nNonZero]
1047       
1048        if (columns) free(columns);                   //delete[] columns;
1049        columns=(int*)malloc(a.nNonZero*sizeof(int)); //columns=new int[a.nNonZero]
1050       
1051        if (this->values)  free(this->values);                   //delete[] values;
1052        this->values=(int*)malloc(a.nNonZero*sizeof(int)); //values=new int[a.nNonZero]
1053     
1054        int  const stop         = this->nNonZero;
1055        T*   const this_values  = this->values;     // MTA -to force onto stack
1056        T*   const a_values     = a.values;         // MTA -to force onto stack
1057        int* const this_rows    = this->rows;       // MTA -to force onto stack
1058        int* const this_columns = this->columns;    // MTA -to force onto stack
1059        int* const a_rows       = a.rows;           // MTA -to force onto stack
1060        int* const a_columns    = a.columns;        // MTA -to force onto stack
1061        #pragma mta assert parallel
1062        for (int i = 0; i < stop; i++)  {
1063          this_values[i]  = a_values[i];
1064          this_rows[i]    = a_rows[i];
1065          this_columns[i] = a_columns[i];
1066        }
1067      }
1068      return *this;
1069    }   
1070   
1071    SparseMatrixCOO<T>& operator= (const SparseMatrixCSR<T>& a) 
1072    {
1073        printf("NOT IMPLEMENTED: Converting from CSR to COO\n");
1074        return *this;
1075    }
1076     
1077    SparseMatrixCOO<T>& operator= (const SparseMatrixCSC<T>& a) 
1078    {
1079        printf("NOT IMPLEMENTED: Converting from CSC to COO\n");
1080        return *this;
1081    }
1082           
1083    SparseMatrixCOO<T>& operator() (const SparseMatrixCSR<T>&);
1084    SparseMatrixCOO<T>& operator() (const SparseMatrixCSC<T>&);
1085   
1086    DenseVector<T> friend
1087     operator* <T> (const SparseMatrixCOO<T>&, const DenseVector<T>&);
1088   
1089    void MatrixPrint (char const* name) 
1090    {
1091        printf("SparseMatrixCOO Print %s row $d col %d\n",name,this->nRow,this->nCol);
1092    }
1093};
1094
1095
1096
1097template <class T> DenseVector<T>
1098operator* (const SparseMatrixCOO<T>& a, const DenseVector<T>& b)
1099{
1100  if (a.nCol != b.length) 
1101  {
1102      printf("INCOMPATIBLE SparseMatrixCoo * DenseVector multiplication\n");
1103      exit(1);
1104  }
1105   
1106  DenseVector<T> temp(b.length);
1107  int const stop        = temp.length;
1108  T*  const temp_values = temp.values;
1109  T   const zero        = T();
1110  #pragma mta assert parallel
1111  for (int i = 0; i < stop; i++)
1112    temp_values[i] = zero;
1113 
1114  int  const end       = a.nNonZero;
1115  T*   const a_values  = a.values;        // MTA -to force onto stack
1116  T*   const b_values  = b.values;        // MTA -to force onto stack
1117  int* const a_rows    = a.rows;          // MTA -to force onto stack
1118  int* const b_rows    = b.rows;          // MTA -to force onto stack
1119  int* const a_columns = a.columns;       // MTA -to force onto stack
1120  #pragma mta assert parallel
1121  for (int i = 0; i < end; i++)
1122    mt_inc (temp_values[a_rows[i]], a_values[i] * b_values[a_columns[i]]);
1123  return temp;
1124}
1125
1126
1127/* biconjugate gradient solver derived from linbcg.c in "Numerical Recipes */
1128/* in C", second edition, Press, Vetterling, Teukolsky, Flannery, pp 86-89 */
1129template <class T> DenseVector<T>&
1130linbcg (const SparseMatrixCSR<T>& A,
1131  DenseVector<T>& x,
1132  const DenseVector<T>& b,
1133  int const itermax,
1134  T& err,
1135  T const tol)
1136{
1137  int const length = A.MatrixRows();
1138  double const bnorm = b.norm2();
1139  double ak, akden, bk, bknum, bkden;
1140
1141  DenseVector<T> p(length), pp(length);
1142  DenseVector<T> r(length), rr(length);
1143  DenseVector<T> z(length), zz(length);
1144  DenseVector<T> d = diagonal(A);
1145 
1146  r  = b - (A * x);
1147  z  = d.asolve (r);
1148  rr = r;
1149 
1150  for (int iter = 0; iter < itermax; iter++)  {
1151    zz = d.asolve (rr);
1152
1153    bknum = z * rr;
1154   
1155    if (iter == 0)  {
1156       p  = z;
1157       pp = zz;       
1158       }
1159    else  {
1160       bk = bknum / bkden;
1161       p  = (bk * p)  + z;
1162       pp = (bk * pp) + zz;
1163       }
1164    bkden =  bknum;
1165
1166    z     = A * p;
1167    akden = z * pp;
1168    ak    = bknum / akden;
1169
1170    zz = Transpose_SMVm(A, pp);
1171
1172    x  += (ak * p);
1173    r  -= (ak * z);
1174    rr -= (ak * zz);
1175
1176    z = d.asolve(r);
1177    err = r.norm2() / bnorm;
1178    if (err < tol)
1179      break;
1180    }
1181  return x;
1182}
1183#endif
1184
1185// EOF
Note: See TracBrowser for help on using the repository browser.