source: trunk/LinAlgImpl/Serial/IpGenTMatrix.cpp @ 2

Last change on this file since 2 was 2, checked in by andreasw, 15 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.5 KB
Line 
1// Copyright (C) 2004, International BusinDess Machines and others.
2// All Rights Reserved.
3// This code is published under the Common Public License.
4//
5// $Id: IpGenTMatrix.cpp 2 2004-10-21 01:03:09Z andreasw $
6//
7// Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
8
9#include "IpGenTMatrix.hpp"
10#include "IpDenseVector.hpp"
11#include "IpBlas.hpp"
12
13namespace Ipopt
14{
15
16  GenTMatrix::GenTMatrix(const GenTMatrixSpace* owner_space)
17      :
18      Matrix(owner_space),
19      owner_space_(owner_space),
20      values_(NULL),
21      initialized_(false)
22  {
23    values_ = owner_space_->AllocateInternalStorage();
24
25    if (Nonzeros() == 0) {
26      initialized_ = true; // I guess ?!? what does this mean ?!?
27    }
28  }
29
30  GenTMatrix::~GenTMatrix()
31  {
32    owner_space_->FreeInternalStorage(values_);
33  }
34
35  void GenTMatrix::SetValues(const Number* Values)
36  {
37    IpBlasDcopy(Nonzeros(), Values, 1, values_, 1);
38    initialized_ = true;
39    ObjectChanged();
40  }
41
42  void GenTMatrix::MultVectorImpl(Number alpha, const Vector &x, Number beta,
43                                  Vector &y) const
44  {
45    //  A few sanity checks
46    DBG_ASSERT(NCols()==x.Dim());
47    DBG_ASSERT(NRows()==y.Dim());
48
49    // Take care of the y part of the addition
50    DBG_ASSERT(initialized_);
51    if( beta!=0.0 ) {
52      y.Scal(beta);
53    }
54    else {
55      y.Set(0.0);  // In case y hasn't been initialized yet
56    }
57
58    // See if we can understand the data
59    const DenseVector* dense_x = dynamic_cast<const DenseVector*>(&x);
60    DBG_ASSERT(dense_x); /* ToDo: Implement others */
61    DenseVector* dense_y = dynamic_cast<DenseVector*>(&y);
62    DBG_ASSERT(dense_y); /* ToDo: Implement others */
63
64    if (dense_x && dense_y) {
65      const Index*  irows=Irows();
66      const Index*  jcols=Jcols();
67      const Number* val=values_;
68      const Number* xvals=dense_x->Values();
69      Number* yvals=dense_y->Values();
70      for(Index i=0; i<Nonzeros(); i++) {
71        yvals[*irows-1] += alpha* (*val) * xvals[*jcols-1];
72        val++;
73        irows++;
74        jcols++;
75      }
76    }
77  }
78
79  void GenTMatrix::TransMultVectorImpl(Number alpha, const Vector &x, Number beta,
80                                       Vector &y) const
81  {
82    //  A few sanity checks
83    DBG_ASSERT(NCols()==y.Dim());
84    DBG_ASSERT(NRows()==x.Dim());
85
86    // Take care of the y part of the addition
87    DBG_ASSERT(initialized_);
88    if( beta!=0.0 ) {
89      y.Scal(beta);
90    }
91    else {
92      y.Set(0.0);  // In case y hasn't been initialized yet
93    }
94
95    // See if we can understand the data
96    const DenseVector* dense_x = dynamic_cast<const DenseVector*>(&x);
97    DBG_ASSERT(dense_x); /* ToDo: Implement others */
98    DenseVector* dense_y = dynamic_cast<DenseVector*>(&y);
99    DBG_ASSERT(dense_y); /* ToDo: Implement others */
100
101    if (dense_x && dense_y) {
102      const Index*  irows=Irows();
103      const Index*  jcols=Jcols();
104      const Number* val=values_;
105      const Number* xvals=dense_x->Values();
106      Number* yvals=dense_y->Values();
107      for(Index i=0; i<Nonzeros(); i++) {
108        yvals[*jcols-1] += alpha* (*val) * xvals[*irows-1];
109        val++;
110        irows++;
111        jcols++;
112      }
113    }
114  }
115
116  void GenTMatrix::PrintImpl(FILE* fp, std::string name, Index indent, std::string prefix) const
117  {
118    fprintf(fp, "\n");
119    for (Index ind=0; ind<indent; ind++) {
120      fprintf(fp, " ");
121    }
122    fprintf(fp, "%sGenTMatrix \"%s\" with %d nonzero elements:\n",
123            prefix.c_str(), name.c_str(), Nonzeros());
124    if (initialized_) {
125      for (Index i=0; i<Nonzeros(); i++) {
126        for (Index ind=0; ind<indent; ind++) {
127          fprintf(fp, " ");
128        }
129        fprintf(fp, "%s%s[%5d,%5d]=%23.16e  (%d)\n", prefix.c_str(), name.c_str(), Irows()[i],
130                Jcols()[i], values_[i], i);
131      }
132    }
133    else {
134      for (Index ind=0; ind<indent; ind++) {
135        fprintf(fp, " ");
136      }
137      fprintf(fp, "%sUninitialized! ", prefix.c_str());
138    }
139  }
140
141  GenTMatrixSpace::GenTMatrixSpace(Index nRows, Index nCols,
142                                   Index nonZeros,
143                                   const Index* iRows, const Index* jCols)
144      :
145      MatrixSpace(nRows, nCols),
146      nonZeros_(nonZeros),
147      iRows_(NULL),
148      jCols_(NULL)
149  {
150    iRows_ = new Index[nonZeros];
151    jCols_ = new Index[nonZeros];
152    for (Index i=0; i<nonZeros; i++) {
153      iRows_[i] = iRows[i];
154      jCols_[i] = jCols[i];
155    }
156  }
157
158  Number* GenTMatrixSpace::AllocateInternalStorage() const
159  {
160    return new Number[Nonzeros()];
161  }
162
163  void GenTMatrixSpace::FreeInternalStorage(Number* values) const
164  {
165    delete [] values;
166  }
167
168
169} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.