ARTS  2.2.66
matpackII.cc
Go to the documentation of this file.
1 /* Copyright (C) 2001-2012
2  Stefan Buehler <sbuehler@ltu.se>
3  Mattias Ekstroem <ekstrom@rss.chalmers.se>
4 
5  This program is free software; you can redistribute it and/or modify it
6  under the terms of the GNU General Public License as published by the
7  Free Software Foundation; either version 2, or (at your option) any
8  later version.
9 
10  This program is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU General Public License for more details.
14 
15  You should have received a copy of the GNU General Public License
16  along with this program; if not, write to the Free Software
17  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
18  USA. */
19 
40 // #include <vector>
41 #include <algorithm>
42 #include <set>
43 #include <iostream> // For debugging.
44 #include <cmath>
45 #include <iterator>
46 #include "matpackII.h"
47 
48 using std::vector;
49 using std::setw;
50 
51 
52 // Simple member Functions
53 // ----------------
54 
57 {
58  return mrr;
59 }
60 
63 {
64  return mcr;
65 }
66 
69 {
70  return *(mcolptr.end()-1);
71 }
72 
73 // Index Operators
74 // ---------------
75 
77 
94 {
95  // Check if indices are valid:
96  assert( 0<=r );
97  assert( 0<=c );
98  assert( r<mrr );
99  assert( c<mcr );
100 
101  // Get index of first data element of this column:
102  Index i = mcolptr[c];
103 
104  // Get index of first data element of next column. (This always works,
105  // because mcolptr has one extra element pointing behind the last
106  // column.)
107  const Index end = mcolptr[c+1];
108 
109  // See if we find an element with the right row index in this
110  // range. We assume that the elements are sorted by ascending row
111  // index. We use the index i as the loop counter, which we have
112  // initialized above.
113  for ( ; i<end; ++i )
114  {
115  Index rowi = mrowind[i];
116  if ( r < rowi )
117  break;
118  if ( r == rowi )
119  return mdata[i];
120  }
121 
122  // If we are here, then the requested data element does not yet
123  // exist.
124 
125  // We have to adjust the array of column pointers. The values
126  // in all columns above the current one have to be increased by 1.
127  for ( vector<Index>::iterator j = mcolptr.begin() + c + 1;
128  j < mcolptr.end();
129  ++j )
130  ++(*j);
131 
132  // We have to insert the new element in *mrowind and *mdata just one
133  // position before the index i. We can use insert to achieve
134  // this. Because they return an iterator to the newly inserted
135  // element, we can return directly from the second call.
136  mrowind.insert( mrowind.begin()+i, r );
137  return *( mdata.insert( mdata.begin()+i, 0 ) );
138 }
139 
141 
153 { return this->ro(r, c); }
154 
156 
171 {
172  // Check if indices are valid:
173  assert( 0<=r );
174  assert( 0<=c );
175  assert( r<mrr );
176  assert( c<mcr );
177 
178  // Get index of first data element of this column:
179  Index begin = mcolptr[c];
180 
181  // Get index of first data element of next column. (This always works,
182  // because mcolptr has one extra element pointing behind the last
183  // column.)
184  const Index end = mcolptr[c+1];
185 
186  // See if we find an element with the right row index in this
187  // range. We assume that the elements are sorted by ascending row
188  // index.
189  for ( Index i=begin; i<end; ++i )
190  {
191  Index rowi = mrowind[i];
192  if ( r < rowi )
193  return 0;
194  if ( r == rowi )
195  return mdata[i];
196  }
197  return 0;
198 }
199 
200 
201 // Constructors
202 // ------------
203 
205 
210  mdata(),
211  mrowind(),
212  mcolptr(),
213  mrr(0),
214  mcr(0)
215 {
216  // Nothing to do here
217 }
218 
220 
239  mdata(),
240  mrowind(),
241  mcolptr(c+1,0),
242  mrr(r),
243  mcr(c)
244 {
245  // Nothing to do here.
246 }
247 
248 
250 
265 {
266  // Check if the row index and the Vector length are valid
267  assert( 0<=r );
268  assert( r<mrr );
269  assert( v.nelem()==mcr );
270 
271  // Calculate number of non-zero elements in Vector v.
272  Index vnnz=0;
273  for (Index i=0; i<v.nelem(); i++)
274  {
275  if (v[i]!=0)
276  {
277  vnnz++;
278  }
279  }
280 
281  // Count number of already existing elements in this row. Create
282  // reference mrowind and mdata vector that copies of the real mrowind and
283  // mdata. Resize the real mrowind and mdata to the correct output size.
284  Index rnnz = vnnz - count(mrowind.begin(),mrowind.end(),r);
285 
286  vector<Index> mrowind_ref(mrowind.size());
287  copy(mrowind.begin(), mrowind.end(), mrowind_ref.begin());
288 
289  vector<Numeric> mdata_ref(mdata.size());
290  copy(mdata.begin(), mdata.end(), mdata_ref.begin());
291 
292  mrowind.resize(mrowind_ref.size()+rnnz);
293  mdata.resize(mdata_ref.size()+rnnz);
294 
295  // Create iterators to the output vectors to keep track of current
296  // positions.
297  vector<Index>::iterator mrowind_it = mrowind.begin();
298  vector<Numeric>::iterator mdata_it = mdata.begin();
299 
300  // Create a variable to store the change to mcolptr for each run
301  Index colptr_mod = 0;
302 
303  // Loop through Vector v and insert the non-zero elements.
304  for (Index i=0; i<v.nelem(); i++)
305  {
306  // Get mdata- and mrowind iterators to start and end of this
307  // (the i:th) reference column.
308  vector<Numeric>::iterator dstart = mdata_ref.begin()+mcolptr[i];
309  vector<Numeric>::iterator dend = mdata_ref.begin()+mcolptr[i+1];
310  vector<Index>::iterator rstart = mrowind_ref.begin()+mcolptr[i];
311  vector<Index>::iterator rend = mrowind_ref.begin()+mcolptr[i+1];
312 
313  // Apply mcolptr change now that we have the iterators to the
314  // data and row indices.
315  mcolptr[i] = colptr_mod;
316 
317  if (v[i]!=0)
318  {
319  // Check if r exist within this column, and get iterator for
320  // mdata
321  vector<Index>::iterator rpos = find(rstart,rend,r);
322  vector<Numeric>::iterator dpos = dstart+(rpos-rstart);
323  if (rpos!=rend)
324  {
325  // The index was found, replace the value in mdata with the
326  // value from v.
327  *dpos = v[i];
328 
329  // Copy this column to the ouput vectors.
330  copy(rstart,rend,mrowind_it);
331  copy(dstart,dend,mdata_it);
332 
333  // Adjust the position iterators accordingly.
334  mrowind_it += rend-rstart;
335  mdata_it += dend-dstart;
336 
337  // Set the mcolptr step, for next loop
338  colptr_mod = mcolptr[i]+(rend-rstart);
339  }
340  else
341  {
342  // The row index was not found, look for the first index
343  // greater than r
344  rpos = find_if(rstart,rend,bind2nd(std::greater<Index>(),r));
345  dpos = dstart+(rpos-rstart);
346 
347  // Copy the first part of the column to the output vector.
348  copy(rstart,rpos,mrowind_it);
349  copy(dstart,dpos,mdata_it);
350 
351  // Make sure mrowind_it and mdata_it points at the first
352  // 'empty' position.
353  mrowind_it += rpos-rstart;
354  mdata_it += dpos-dstart;
355 
356  // Insert the new value from v in mdata and the row index r in
357  // mrowind.
358  *mrowind_it = r;
359  *mdata_it = v[i];
360 
361  // Again, make sure mrowind_it and mdata_it points at the
362  // first 'empty' position.
363  mrowind_it++;
364  mdata_it++;
365 
366  // Copy the rest of this column to the output vectors.
367  copy(rpos,rend,mrowind_it);
368  copy(dpos,dend,mdata_it);
369 
370  // Adjust the iterators a last time
371  mrowind_it += rend-rpos;
372  mdata_it += dend-dpos;
373 
374  // Set the mcolptr step, for next loop
375  colptr_mod = mcolptr[i]+(rend-rstart+1);
376  }
377  }
378  else
379  {
380  // Check if r exist within this column, and get iterator for
381  // mdata
382  vector<Index>::iterator rpos = find(rstart,rend,r);
383  vector<Numeric>::iterator dpos = dstart+(rpos-rstart);
384  if (rpos!=rend)
385  {
386  // The index was found and we use remove_copy to copy the rest
387  // of the column to mrowind_new and mdata_new.
388  remove_copy(rstart,rend,mrowind_it, *rpos);
389  remove_copy(dstart,dend,mdata_it, *dpos);
390 
391  // Increase the mrowind_it and mdata_it
392  mrowind_it += rend-rstart-1;
393  mdata_it += dend-dstart-1;
394 
395  // Set the mcolptr step, for next loop
396  colptr_mod = mcolptr[i]+(rend-rstart-1);
397  }
398  else
399  {
400  // The row index was not found, all we need is to copy this
401  // columns to the output mrowind and mdata vectors.
402  copy(rstart,rend,mrowind_it);
403  copy(dstart,dend,mdata_it);
404 
405  // Adjust the mrowind_it and mdata_it iterators
406  mrowind_it += rend-rstart;
407  mdata_it += dend-dstart;
408 
409  // Set the mcolptr step, for next loop
410  colptr_mod = mcolptr[i]+(rend-rstart);
411  }
412  }
413  }
414  // Apply mcolptr change for the one extra mcolptr element.
415  *(mcolptr.end()-1) = colptr_mod;
416 }
417 
419 
431 {
432  assert( 0<=r );
433  assert( 0<=c );
434 
435  // First get number of ones in the identity matrix
436  Index n = std::min(r,c);
437 
438  // Remake and assign values to vectors
439  mcolptr = vector<Index>( c+1, n);
440  mrowind = vector<Index>(n);
441  mdata = vector<Numeric>(n,1.0);
442 
443  // Loop over number of ones and assign values [0:n-1] to mcolptr and
444  // mrowind
445  vector<Index>::iterator rit = mrowind.begin();
446  vector<Index>::iterator cit = mcolptr.begin();
447  for( Index i=0; i<n; i++ ) {
448  *rit++ = i;
449  *cit++ = i;
450  }
451 
452  mrr = r;
453  mcr = c;
454 }
455 
457 
470 {
471  assert( 0<=r );
472  assert( 0<=c );
473  if ( mrr!=r || mcr!=c )
474  {
475  mdata.resize(0);
476  mrowind.resize(0);
477  mcolptr = vector<Index>(c+1,0);
478 
479  mrr = r;
480  mcr = c;
481  }
482 }
483 
484 
486 
495 std::ostream& operator<<(std::ostream& os, const Sparse& v)
496 {
497  for (size_t c=0; c<v.mcolptr.size()-1; ++c)
498  {
499  // Get index of first data element of this column:
500  Index begin = v.mcolptr[c];
501 
502  // Get index of first data element of next column. (This always works,
503  // because mcolptr has one extra element pointing behind the last
504  // column.)
505  const Index end = v.mcolptr[c+1];
506 
507  // Loop through the elements in this column:
508  for ( Index i=begin; i<end; ++i )
509  {
510  // Get row index:
511  Index r = v.mrowind[i];
512 
513  // Output everything:
514  os << setw(3) << r << " "
515  << setw(3) << c << " "
516  << setw(3) << v.mdata[i] << "\n";
517  }
518  }
519 
520  return os;
521 }
522 
523 // General matrix functions
524 
526 
537 void abs( Sparse& A,
538  const Sparse& B )
539 {
540  // Check dimensions
541  assert( A.nrows() == B.nrows() );
542  assert( A.ncols() == B.ncols() );
543 
544  // Here we allocate memory for the A.mdata vector so that it matches the
545  // input matrix, and then store the absolute values of B.mdata in it.
546  A.mdata.resize( B.mdata.size() );
547  Index end = B.mdata.size();
548  for (Index i=0; i<end; i++)
549  {
550  A.mdata[i] = fabs(B.mdata[i]);
551  }
552 
553  // The column pointer and row index vectors are copies of the input
554  A.mcolptr.resize( B.mcolptr.size() );
555  copy(B.mcolptr.begin(), B.mcolptr.end(), A.mcolptr.begin());
556  A.mrowind.resize( B.mrowind.size() );
557  copy(B.mrowind.begin(), B.mrowind.end(), A.mrowind.begin());
558 }
559 
560 
562 
579 void mult( VectorView y,
580  const Sparse& M,
581  ConstVectorView x )
582 {
583  // Check dimensions:
584  assert( y.nelem() == M.nrows() );
585  assert( M.ncols() == x.nelem() );
586 
587  // Initialize y to all zeros:
588  y = 0.0;
589 
590  // Looping through every element of M, multiplying it with the
591  // appropriate elements of x.
592  for (size_t c=0; c<M.mcolptr.size()-1; ++c)
593  {
594  // Get index of first data element of this column:
595  Index begin = M.mcolptr[c];
596 
597  // Get index of first data element of next column. (This always works,
598  // because mcolptr has one extra element pointing behind the last
599  // column.)
600  const Index end = M.mcolptr[c+1];
601 
602  // Loop through the elements in this column:
603  for ( Index i=begin; i<end; ++i )
604  {
605  // Get row index:
606  Index r = M.mrowind[i];
607 
608  // Compute this element:
609  y[r] += M.mdata[i] * x[c];
610  }
611  }
612 }
613 
614 
616 
633 void mult( MatrixView A,
634  const Sparse& B,
636 {
637  // Check dimensions:
638  assert( A.nrows() == B.nrows() );
639  assert( A.ncols() == C.ncols() );
640  assert( B.ncols() == C.nrows() );
641 
642  // Set all elements of A to zero:
643  A = 0.0;
644 
645  // Loop through the elements of B:
646  for (size_t c=0; c<B.mcolptr.size()-1; ++c)
647  {
648  // Get index of first data element of this column:
649  Index begin = B.mcolptr[c];
650 
651  // Get index of first data element of next column. (This always works,
652  // because mcolptr has one extra element pointing behind the last
653  // column.)
654  const Index end = B.mcolptr[c+1];
655 
656  // Loop through the elements in this column:
657  for ( Index i=begin; i<end; ++i )
658  {
659  // Get row index:
660  Index r = B.mrowind[i];
661 
662  // Multiply this element with the corresponding row of C and
663  // add the product to the right row of A
664  for (Index j=0; j<C.ncols(); j++)
665  {
666  /* Conceptually:
667  A(r,j) += B.ro(r,c) * C(c,j);
668  But we don't need to use the index operator, because
669  we have the right element right here: */
670  A(r,j) += B.mdata[i] * C(c,j);
671  }
672  }
673  }
674 }
675 
676 
678 
689 void transpose( Sparse& A,
690  const Sparse& B )
691 {
692  // Check dimensions
693  assert( A.nrows() == B.ncols() );
694  assert( A.ncols() == B.nrows() );
695  if ( B.mdata.size() == 0) return;
696 
697  // Here we allocate memory for the A.mdata and A.mrowind vectors so that
698  // it matches the input matrix. (NB: that mdata and mrowind already has
699  // the same size.)
700  A.mdata.resize( B.mdata.size() );
701  A.mrowind.resize( B.mrowind.size() );
702 
703  // Find the minimum and maximum existing row number in B.
704  // (This maybe unnecessary in our case since we mostly treat diagonally
705  // banded matrices where the minimum and maximum existing row numbers
706  // coincide with the border numbers of the matrix.)
707  vector<Index>::const_iterator startrow, stoprow;
708  startrow = min_element( B.mrowind.begin(), B.mrowind.end() );
709  stoprow = max_element( B.mrowind.begin(), B.mrowind.end() );
710 
711  // To create the A.mcolptr vector we loop through the existing row
712  // numbers of B and count how many there are in B.mrowind.
713  // First we make sure it is initialized to all zeros.
714  A.mcolptr.assign( A.mcr+1, 0 );
715  for (Index i=*startrow; i<=*stoprow; i++)
716  {
717  Index n = count( B.mrowind.begin(), B.mrowind.end(), i);
718 
719  // The column pointer for the column above the current one in A
720  // are then set to the current one plus this amount.
721  A.mcolptr[i+1] = A.mcolptr[i] + n;
722  }
723 
724  // Next we loop through the columns of A and search for the corresponding
725  // row index within each column of B (i.e. two loops). The elements are
726  // then stored in A. We know that for every column in B we will have
727  // the corresponding rowindex in A.
728  vector<Numeric>::iterator dataA_it = A.mdata.begin();
729  vector<Index>::iterator rowindA_it = A.mrowind.begin();
730 
731  // Looping over columns in A to keep track of what row number we are
732  // looking for.
733  for (size_t c=0; c<A.mcolptr.size(); ++c)
734  {
735  // Looping over columns in B to get the elements in the right order,
736  // this will turn into row index in A.
737  for (size_t i=0; i<B.mcolptr.size()-1; i++)
738  {
739  // Since only one element can occupy one row in a column, we
740  // only need to call find() once per column in B.
741  vector<Index>::const_iterator elem_it =
742  find(B.mrowind.begin()+*(B.mcolptr.begin()+i),
743  B.mrowind.begin()+*(B.mcolptr.begin()+i+1), Index (c));
744 
745  // If we found the element, store it in A.mrowind and A.mdata
746  if (elem_it != B.mrowind.begin()+*(B.mcolptr.begin()+i+1))
747  {
748  *rowindA_it = i;
749  rowindA_it++;
750 
751  // To get the corresponding element in B.mdata we subtract
752  // the initial value of the row index iterator from elem_it
753  // and add the difference to a mdata iterator.
754  Index diff = elem_it - B.mrowind.begin();
755  vector<Numeric>::const_iterator elemdata_it=B.mdata.begin() + diff;
756  *dataA_it = *elemdata_it;
757  dataA_it++;
758  }
759  }
760  }
761 }
762 
763 
765 
780 void mult( Sparse& A,
781  const Sparse& B,
782  const Sparse& C )
783 {
784  // Check dimensions and make sure that A is empty
785  assert( A.nrows() == B.nrows() );
786  assert( A.ncols() == C.ncols() );
787  assert( B.ncols() == C.nrows() );
788  A.mcolptr.assign( A.mcr+1, 0);
789  A.mrowind.clear();
790  A.mdata.clear();
791 
792  // Transpose B to simplify multiplication algorithm, after transposing we
793  // can extract columns form the two matrices and multiply them, (which is
794  // easier than extacting rows.)
795  Sparse Bt(B.ncols(), B.nrows());
796  // cout << "B.nrows, B.ncols = " << B.nrows() << ", " << B.ncols() << "\n";
797  // cout << "B = \n" << B << "\n";
798  transpose(Bt,B);
799 
800  // By looping over columns in C and multiply them with every column in Bt
801  // (instead of the conventional loooping over Bt), we get the output
802  // elements in the right order for storing them.
803  for (size_t c=0; c<C.mcolptr.size()-1; ++c)
804  {
805  // Get row indices of this column
806  //Index beginC = (*C.mcolptr)[c];
807  //Index endC = (*C.mcolptr)[c+1];
808 
809  for (size_t b=0; b<Bt.mcolptr.size()-1; ++b)
810  {
811  // Get the intersection between the elements in the two columns and
812  // and store them in a temporary vector
813  std::set<Index> colintersec;
814  set_intersection(C.mrowind.begin()+*(C.mcolptr.begin()+c),
815  C.mrowind.begin()+*(C.mcolptr.begin()+c+1),
816  Bt.mrowind.begin()+*(Bt.mcolptr.begin()+b),
817  Bt.mrowind.begin()+*(Bt.mcolptr.begin()+b+1),
818  inserter(colintersec, colintersec.begin()));
819 
820  // If we got an intersection, loop through it and multiply the
821  // element pairs from C and Bt and store result in A
822  if (!colintersec.empty())
823  {
824  Numeric tempA = 0.0;
825  for (std::set<Index>::iterator i=colintersec.begin();
826  i!=colintersec.end(); ++i)
827  {
828  // To get iterators to the data elements in Bt and C, we
829  // subtract the iterator that points to the start of the
830  // mrowind vector from the iterators that points to the
831  // values in colintersec. We then get the number of steps
832  // from the start of the both vectors mrowind and mdata to
833  // the intersecting values and can add them to the
834  // iterator that points to the beginning of mdata.
835  vector<Index>::const_iterator rowindBt_it =
836  find(Bt.mrowind.begin()+*(Bt.mcolptr.begin()+b),
837  Bt.mrowind.begin()+*(Bt.mcolptr.begin()+b+1), *i);
838  vector<Index>::const_iterator rowindC_it =
839  find(C.mrowind.begin()+*(C.mcolptr.begin()+c),
840  C.mrowind.begin()+*(C.mcolptr.begin()+c+1), *i);
841 
842  vector<Numeric>::const_iterator dataBt_it =
843  Bt.mdata.begin()+(rowindBt_it-Bt.mrowind.begin());
844  vector<Numeric>::const_iterator dataC_it =
845  C.mdata.begin()+(rowindC_it-C.mrowind.begin());
846 
847  tempA += *dataBt_it * *dataC_it;
848  }
849  A.rw(b,c) = tempA;
850  }
851  }
852  }
853 }
854 
855 
857 
871 void add( Sparse& A,
872  const Sparse& B,
873  const Sparse& C )
874 {
875  // Check dimensions
876  assert( B.ncols() == C.ncols() );
877  assert( B.nrows() == C.nrows() );
878 
879  // We copy the smaller matrix of B or C to A. This way we can
880  // loop over the matrix with the fewer number of elements to perform
881  // the actual addition later.
882  const Sparse* D;
883  if (B.data().size() < C.data().size())
884  {
885  A=C;
886  D=&B;
887  }
888  else
889  {
890  A=B;
891  D=&C;
892  }
893 
894  for (size_t c = 0; c < D->mcolptr.size()-1; ++c)
895  {
896  // Loop through the elements in this column:
897  for (Index i = D->mcolptr[c]; i < D->mcolptr[c+1]; ++i)
898  {
899  A.rw(D->mrowind[i], c) += D->mdata[i];
900  }
901  }
902 }
903 
904 
906 
920 void sub( Sparse& A,
921  const Sparse& B,
922  const Sparse& C )
923 {
924  // Check dimensions
925  assert( B.ncols() == C.ncols() );
926  assert( B.nrows() == C.nrows() );
927 
928  A=B;
929 
930  for (size_t c = 0; c < C.mcolptr.size()-1; ++c)
931  {
932  // Loop through the elements in this column:
933  for (Index i = C.mcolptr[c]; i < C.mcolptr[c+1]; ++i)
934  {
935  A.rw(C.mrowind[i], c) -= C.mdata[i];
936  }
937  }
938 }
939 
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
The VectorView class.
Definition: matpackI.h:372
Index nnz() const
Returns the number of nonzero elements.
Definition: matpackII.cc:68
Numeric operator()(Index r, Index c) const
Plain index operator.
Definition: matpackII.cc:152
#define C(a, b)
Definition: Faddeeva.cc:254
#define M
Definition: rng.cc:196
Index mrr
Number of rows in the sparse matrix.
Definition: matpackII.h:102
The Vector class.
Definition: matpackI.h:556
The MatrixView class.
Definition: matpackI.h:679
void copy(ConstIterator1D origin, const ConstIterator1D &end, Iterator1D target)
Copy data between begin and end to target.
Definition: matpackI.cc:612
The Sparse class.
Definition: matpackII.h:55
const std::vector< Numeric > & data() const
Definition: matpackII.h:75
Index mcr
Number of rows in the sparse matrix.
Definition: matpackII.h:104
Index ncols() const
Returns the number of columns.
Definition: matpackII.cc:62
Numeric & rw(Index r, Index c)
Read and write index operator.
Definition: matpackII.cc:93
void make_I(Index r, Index c)
Make Identity matrix.
Definition: matpackII.cc:430
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:838
Sparse()
Default constructor.
Definition: matpackII.cc:209
Header file for sparse matrices.
friend void transpose(Sparse &A, const Sparse &B)
Transpose of sparse matrix.
Definition: matpackII.cc:689
friend void abs(Sparse &A, const Sparse &B)
Absolute value of sparse matrix elements.
Definition: matpackII.cc:537
Index nrows() const
Returns the number of rows.
Definition: matpackII.cc:56
friend void mult(VectorView y, const Sparse &M, ConstVectorView x)
Sparse matrix - Vector multiplication.
Definition: matpackII.cc:579
std::vector< Index > mcolptr
Pointers to first data element for each column.
Definition: matpackII.h:100
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
friend void sub(Sparse &A, const Sparse &B, const Sparse &C)
Sparse - Sparse subtraction.
Definition: matpackII.cc:920
friend std::ostream & operator<<(std::ostream &os, const Sparse &v)
Output operator for Sparse.
Definition: matpackII.cc:495
#define min(a, b)
Definition: continua.cc:20460
std::vector< Index > mrowind
Row indices.
Definition: matpackII.h:98
Numeric ro(Index r, Index c) const
Read only index operator.
Definition: matpackII.cc:170
A constant view of a Vector.
Definition: matpackI.h:292
A constant view of a Matrix.
Definition: matpackI.h:596
void resize(Index r, Index c)
Resize function.
Definition: matpackII.cc:469
std::vector< Numeric > mdata
The actual data values.
Definition: matpackII.h:96
friend void add(Sparse &A, const Sparse &B, const Sparse &C)
Sparse - Sparse addition.
Definition: matpackII.cc:871
void insert_row(Index r, Vector v)
Insert row function.
Definition: matpackII.cc:264
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:832