source: sasview/realspace_modeling/iqPy/libiqPy/tnt/tnt_cmat.h @ a778da5

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since a778da5 was a778da5, checked in by Gervaise Alina <gervyh@…>, 13 years ago

rename folder to allow lower case only and underscore

  • Property mode set to 100644
File size: 10.5 KB
Line 
1/*
2*
3* Template Numerical Toolkit (TNT)
4*
5* Mathematical and Computational Sciences Division
6* National Institute of Technology,
7* Gaithersburg, MD USA
8*
9*
10* This software was developed at the National Institute of Standards and
11* Technology (NIST) by employees of the Federal Government in the course
12* of their official duties. Pursuant to title 17 Section 105 of the
13* United States Code, this software is not subject to copyright protection
14* and is in the public domain. NIST assumes no responsibility whatsoever for
15* its use by other parties, and makes no guarantees, expressed or implied,
16* about its quality, reliability, or any other characteristic.
17*
18*/
19
20
21// C compatible matrix: row-oriented, 0-based [i][j] and 1-based (i,j) indexing
22//
23
24#ifndef TNT_CMAT_H
25#define TNT_CMAT_H
26
27#include "tnt_subscript.h"
28#include "tnt_vec.h"
29#include <cstdlib>
30#include <cassert>
31#include <iostream>
32#include <sstream>
33
34namespace TNT
35{
36
37
38template <class T>
39class Matrix
40{
41
42
43  public:
44
45    typedef Subscript   size_type;
46    typedef         T   value_type;
47    typedef         T   element_type;
48    typedef         T*  pointer;
49    typedef         T*  iterator;
50    typedef         T&  reference;
51    typedef const   T*  const_iterator;
52    typedef const   T&  const_reference;
53
54    Subscript lbound() const { return 1;}
55 
56  protected:
57    Subscript m_;
58    Subscript n_;
59    Subscript mn_;      // total size
60    T* v_;                 
61    T** row_;           
62    T* vm1_ ;       // these point to the same data, but are 1-based
63    T** rowm1_;
64
65    // internal helper function to create the array
66    // of row pointers
67
68    void initialize(Subscript M, Subscript N)
69    {
70        mn_ = M*N;
71        m_ = M;
72        n_ = N;
73
74        v_ = new T[mn_]; 
75        row_ = new T*[M];
76        rowm1_ = new T*[M];
77
78        assert(v_  != NULL);
79        assert(row_  != NULL);
80        assert(rowm1_ != NULL);
81
82        T* p = v_;             
83        vm1_ = v_ - 1;
84        for (Subscript i=0; i<M; i++)
85        {
86            row_[i] = p;
87            rowm1_[i] = p-1;
88            p += N ;
89           
90        }
91
92        rowm1_ -- ;     // compensate for 1-based offset
93    }
94   
95    void copy(const T*  v)
96    {
97        Subscript N = m_ * n_;
98        Subscript i;
99
100#ifdef TNT_UNROLL_LOOPS
101        Subscript Nmod4 = N & 3;
102        Subscript N4 = N - Nmod4;
103
104        for (i=0; i<N4; i+=4)
105        {
106            v_[i] = v[i];
107            v_[i+1] = v[i+1];
108            v_[i+2] = v[i+2];
109            v_[i+3] = v[i+3];
110        }
111
112        for (i=N4; i< N; i++)
113            v_[i] = v[i];
114#else
115
116        for (i=0; i< N; i++)
117            v_[i] = v[i];
118#endif     
119    }
120
121    void set(const T& val)
122    {
123        Subscript N = m_ * n_;
124        Subscript i;
125
126#ifdef TNT_UNROLL_LOOPS
127        Subscript Nmod4 = N & 3;
128        Subscript N4 = N - Nmod4;
129
130        for (i=0; i<N4; i+=4)
131        {
132            v_[i] = val;
133            v_[i+1] = val;
134            v_[i+2] = val;
135            v_[i+3] = val; 
136        }
137
138        for (i=N4; i< N; i++)
139            v_[i] = val;
140#else
141
142        for (i=0; i< N; i++)
143            v_[i] = val;
144       
145#endif     
146    }
147   
148
149   
150    void destroy()
151    {     
152        /* do nothing, if no memory has been previously allocated */
153        if (v_ == NULL) return ;
154
155        /* if we are here, then matrix was previously allocated */
156        if (v_ != NULL) delete [] (v_);     
157        if (row_ != NULL) delete [] (row_);
158
159        /* return rowm1_ back to original value */
160        rowm1_ ++;
161        if (rowm1_ != NULL ) delete [] (rowm1_);
162    }
163
164
165  public:
166
167    operator T**(){ return  row_; }
168    operator T**() const { return row_; }
169
170
171    Subscript size() const { return mn_; }
172
173    // constructors
174
175    Matrix() : m_(0), n_(0), mn_(0), v_(0), row_(0), vm1_(0), rowm1_(0) {};
176
177    Matrix(const Matrix<T> &A)
178    {
179        initialize(A.m_, A.n_);
180        copy(A.v_);
181    }
182
183    Matrix(Subscript M, Subscript N, const T& value = T())
184    {
185        initialize(M,N);
186        set(value);
187    }
188
189    Matrix(Subscript M, Subscript N, const T* v)
190    {
191        initialize(M,N);
192        copy(v);
193    }
194
195    Matrix(Subscript M, Subscript N, const char *s)
196    {
197        initialize(M,N);
198        //std::istrstream ins(s);
199        std::istringstream ins(s);
200
201        Subscript i, j;
202
203        for (i=0; i<M; i++)
204            for (j=0; j<N; j++)
205                ins >> row_[i][j];
206    }
207
208    // destructor
209    //
210    ~Matrix()
211    {
212        destroy();
213    }
214
215
216    // reallocating
217    //
218    Matrix<T>& newsize(Subscript M, Subscript N)
219    {
220        if (num_rows() == M && num_cols() == N)
221            return *this;
222
223        destroy();
224        initialize(M,N);
225       
226        return *this;
227    }
228
229
230
231
232    // assignments
233    //
234    Matrix<T>& operator=(const Matrix<T> &A)
235    {
236        if (v_ == A.v_)
237            return *this;
238
239        if (m_ == A.m_  && n_ == A.n_)      // no need to re-alloc
240            copy(A.v_);
241
242        else
243        {
244            destroy();
245            initialize(A.m_, A.n_);
246            copy(A.v_);
247        }
248
249        return *this;
250    }
251       
252    Matrix<T>& operator=(const T& scalar)
253    { 
254        set(scalar); 
255        return *this;
256    }
257
258
259    Subscript dim(Subscript d) const 
260    {
261#ifdef TNT_BOUNDS_CHECK
262        assert( d >= 1);
263        assert( d <= 2);
264#endif
265        return (d==1) ? m_ : ((d==2) ? n_ : 0); 
266    }
267
268    Subscript num_rows() const { return m_; }
269    Subscript num_cols() const { return n_; }
270
271
272
273
274    inline T* operator[](Subscript i)
275    {
276#ifdef TNT_BOUNDS_CHECK
277        assert(0<=i);
278        assert(i < m_) ;
279#endif
280        return row_[i];
281    }
282
283    inline const T* operator[](Subscript i) const
284    {
285#ifdef TNT_BOUNDS_CHECK
286        assert(0<=i);
287        assert(i < m_) ;
288#endif
289        return row_[i];
290    }
291
292    inline reference operator()(Subscript i)
293    { 
294#ifdef TNT_BOUNDS_CHECK
295        assert(1<=i);
296        assert(i <= mn_) ;
297#endif
298        return vm1_[i]; 
299    }
300
301    inline const_reference operator()(Subscript i) const
302    { 
303#ifdef TNT_BOUNDS_CHECK
304        assert(1<=i);
305        assert(i <= mn_) ;
306#endif
307        return vm1_[i]; 
308    }
309
310
311
312    inline reference operator()(Subscript i, Subscript j)
313    { 
314#ifdef TNT_BOUNDS_CHECK
315        assert(1<=i);
316        assert(i <= m_) ;
317        assert(1<=j);
318        assert(j <= n_);
319#endif
320        return  rowm1_[i][j]; 
321    }
322
323
324   
325    inline const_reference operator() (Subscript i, Subscript j) const
326    {
327#ifdef TNT_BOUNDS_CHECK
328        assert(1<=i);
329        assert(i <= m_) ;
330        assert(1<=j);
331        assert(j <= n_);
332#endif
333        return rowm1_[i][j]; 
334    }
335
336
337
338
339};
340
341
342/* ***************************  I/O  ********************************/
343
344template <class T>
345std::ostream& operator<<(std::ostream &s, const Matrix<T> &A)
346{
347    Subscript M=A.num_rows();
348    Subscript N=A.num_cols();
349
350    s << M << " " << N << "\n";
351
352    for (Subscript i=0; i<M; i++)
353    {
354        for (Subscript j=0; j<N; j++)
355        {
356            s << A[i][j] << " ";
357        }
358        s << "\n";
359    }
360
361
362    return s;
363}
364
365template <class T>
366std::istream& operator>>(std::istream &s, Matrix<T> &A)
367{
368
369    Subscript M, N;
370
371    s >> M >> N;
372
373    if ( !(M == A.num_rows() && N == A.num_cols() ))
374    {
375        A.newsize(M,N);
376    }
377
378
379    for (Subscript i=0; i<M; i++)
380        for (Subscript j=0; j<N; j++)
381        {
382            s >>  A[i][j];
383        }
384
385
386    return s;
387}
388
389// *******************[ basic matrix algorithms ]***************************
390
391
392template <class T>
393Matrix<T> operator+(const Matrix<T> &A, 
394    const Matrix<T> &B)
395{
396    Subscript M = A.num_rows();
397    Subscript N = A.num_cols();
398
399    assert(M==B.num_rows());
400    assert(N==B.num_cols());
401
402    Matrix<T> tmp(M,N);
403    Subscript i,j;
404
405    for (i=0; i<M; i++)
406        for (j=0; j<N; j++)
407            tmp[i][j] = A[i][j] + B[i][j];
408
409    return tmp;
410}
411
412template <class T>
413Matrix<T> operator-(const Matrix<T> &A, 
414    const Matrix<T> &B)
415{
416    Subscript M = A.num_rows();
417    Subscript N = A.num_cols();
418
419    assert(M==B.num_rows());
420    assert(N==B.num_cols());
421
422    Matrix<T> tmp(M,N);
423    Subscript i,j;
424
425    for (i=0; i<M; i++)
426        for (j=0; j<N; j++)
427            tmp[i][j] = A[i][j] - B[i][j];
428
429    return tmp;
430}
431
432template <class T>
433Matrix<T> mult_element(const Matrix<T> &A, 
434    const Matrix<T> &B)
435{
436    Subscript M = A.num_rows();
437    Subscript N = A.num_cols();
438
439    assert(M==B.num_rows());
440    assert(N==B.num_cols());
441
442    Matrix<T> tmp(M,N);
443    Subscript i,j;
444
445    for (i=0; i<M; i++)
446        for (j=0; j<N; j++)
447            tmp[i][j] = A[i][j] * B[i][j];
448
449    return tmp;
450}
451
452
453template <class T>
454Matrix<T> transpose(const Matrix<T> &A)
455{
456    Subscript M = A.num_rows();
457    Subscript N = A.num_cols();
458
459    Matrix<T> S(N,M);
460    Subscript i, j;
461
462    for (i=0; i<M; i++)
463        for (j=0; j<N; j++)
464            S[j][i] = A[i][j];
465
466    return S;
467}
468
469
470   
471template <class T>
472inline Matrix<T> matmult(const Matrix<T>  &A, 
473    const Matrix<T> &B)
474{
475
476#ifdef TNT_BOUNDS_CHECK
477    assert(A.num_cols() == B.num_rows());
478#endif
479
480    Subscript M = A.num_rows();
481    Subscript N = A.num_cols();
482    Subscript K = B.num_cols();
483
484    Matrix<T> tmp(M,K);
485    T sum;
486
487    for (Subscript i=0; i<M; i++)
488    for (Subscript k=0; k<K; k++)
489    {
490        sum = 0;
491        for (Subscript j=0; j<N; j++)
492            sum = sum +  A[i][j] * B[j][k];
493
494        tmp[i][k] = sum; 
495    }
496
497    return tmp;
498}
499
500template <class T>
501inline Matrix<T> operator*(const Matrix<T>  &A, 
502    const Matrix<T> &B)
503{
504    return matmult(A,B);
505}
506
507template <class T>
508inline int matmult(Matrix<T>& C, const Matrix<T>  &A, 
509    const Matrix<T> &B)
510{
511
512    assert(A.num_cols() == B.num_rows());
513
514    Subscript M = A.num_rows();
515    Subscript N = A.num_cols();
516    Subscript K = B.num_cols();
517
518    C.newsize(M,K);
519
520    T sum;
521
522    const T* row_i;
523    const T* col_k;
524
525    for (Subscript i=0; i<M; i++)
526    for (Subscript k=0; k<K; k++)
527    {
528        row_i  = &(A[i][0]);
529        col_k  = &(B[0][k]);
530        sum = 0;
531        for (Subscript j=0; j<N; j++)
532        {
533            sum  += *row_i * *col_k;
534            row_i++;
535            col_k += K;
536        }
537        C[i][k] = sum; 
538    }
539
540    return 0;
541}
542
543
544template <class T>
545Vector<T> matmult(const Matrix<T>  &A, const Vector<T> &x)
546{
547
548#ifdef TNT_BOUNDS_CHECK
549    assert(A.num_cols() == x.dim());
550#endif
551
552    Subscript M = A.num_rows();
553    Subscript N = A.num_cols();
554
555    Vector<T> tmp(M);
556    T sum;
557
558    for (Subscript i=0; i<M; i++)
559    {
560        sum = 0;
561        const T* rowi = A[i];
562        for (Subscript j=0; j<N; j++)
563            sum = sum +  rowi[j] * x[j];
564
565        tmp[i] = sum; 
566    }
567
568    return tmp;
569}
570
571template <class T>
572inline Vector<T> operator*(const Matrix<T>  &A, const Vector<T> &x)
573{
574    return matmult(A,x);
575}
576
577} // namespace TNT
578
579#endif
580// CMAT_H
Note: See TracBrowser for help on using the repository browser.