1 /*************************************************************************
2 * *
3 * Vega FEM Simulation Library Version 3.1 *
4 * *
5 * "sparseMatrix" library , Copyright (C) 2007 CMU, 2009 MIT, 2016 USC *
6 * All rights reserved. *
7 * *
8 * Code author: Jernej Barbic *
9 * http://www.jernejbarbic.com/code *
10 * *
11 * Research: Jernej Barbic, Fun Shing Sin, Daniel Schroeder, *
12 * Doug L. James, Jovan Popovic *
13 * *
14 * Funding: National Science Foundation, Link Foundation, *
15 * Singapore-MIT GAMBIT Game Lab, *
16 * Zumberge Research and Innovation Fund at USC *
17 * *
18 * This library is free software; you can redistribute it and/or *
19 * modify it under the terms of the BSD-style license that is *
20 * included with this library in the file LICENSE.txt *
21 * *
22 * This library is distributed in the hope that it will be useful, *
23 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
24 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the file *
25 * LICENSE.TXT for more details. *
26 * *
27 *************************************************************************/
28
29 #include <stdlib.h>
30 #include <stdio.h>
31 #include <string.h>
32 #include <math.h>
33 #include <cassert>
34 #include "sparseMatrix.h"
35 using namespace std;
36
SparseMatrixOutline(int numRows_)37 SparseMatrixOutline::SparseMatrixOutline(int numRows_): numRows(numRows_)
38 {
39 Allocate();
40 }
41
~SparseMatrixOutline()42 SparseMatrixOutline::~SparseMatrixOutline()
43 {
44 // deallocate column entries
45 for(int i=0; i<numRows; i++)
46 columnEntries[i].clear();
47 columnEntries.clear();
48 }
49
SparseMatrixOutline(int numRows_,double diagonal)50 SparseMatrixOutline::SparseMatrixOutline(int numRows_, double diagonal): numRows(numRows_)
51 {
52 Allocate();
53
54 pair<int,double> entry;
55
56 for(int i=0; i<numRows; i++)
57 {
58 entry.first = i;
59 entry.second = diagonal;
60 columnEntries[i].insert(entry);
61 }
62 }
63
SparseMatrixOutline(int numRows_,double * diagonal)64 SparseMatrixOutline::SparseMatrixOutline(int numRows_, double * diagonal): numRows(numRows_)
65 {
66 Allocate();
67 pair<int,double> entry;
68 for(int i=0; i<numRows; i++)
69 {
70 entry.first = i;
71 entry.second = diagonal[i];
72 columnEntries[i].insert(entry);
73 }
74 }
75
SparseMatrixOutline(const char * filename,int expand)76 SparseMatrixOutline::SparseMatrixOutline(const char * filename, int expand)
77 {
78 if (expand <= 0)
79 {
80 printf("Error: invalid expand factor %d in SparseMatrixOutline constructor.\n", expand);
81 throw 1;
82 }
83
84 FILE * inputMatrix = fopen(filename,"r");
85 if (!inputMatrix)
86 {
87 printf("Error: couldn't open input sparse matrix file %s.\n", filename);
88 throw 2;
89 }
90
91 // read input size
92 int m1, n1;
93 if (fscanf(inputMatrix, "%d\n%d\n", &m1, &n1) < 2)
94 {
95 printf("Error: could not read sparse matrix dimensions in file %s.\n", filename);
96 throw 3;
97 }
98
99 numRows = expand * m1;
100
101 printf("Loading matrix from %s... Size is %d x %d .\n", filename, numRows, expand * n1);fflush(NULL);
102
103 Allocate();
104
105 char s[4096];
106 while (fgets(s,4096,inputMatrix) != NULL)
107 {
108 int i1,j1;
109 double x;
110 sscanf(s, "%d %d %lf\n", &i1, &j1, &x);
111 for(int e=0; e<expand; e++)
112 AddEntry(expand * i1 + e, expand * j1 + e, x);
113 }
114
115 fclose(inputMatrix);
116 }
117
Allocate()118 void SparseMatrixOutline::Allocate()
119 {
120 // allocate empty datastructure for row entries
121 columnEntries.clear();
122 columnEntries.resize(numRows);
123 }
124
IncreaseNumRows(int numAddedRows)125 void SparseMatrixOutline::IncreaseNumRows(int numAddedRows)
126 {
127 numRows += numAddedRows;
128 columnEntries.resize(numRows);
129 }
130
AddEntry(int i,int j,double value)131 void SparseMatrixOutline::AddEntry(int i, int j, double value)
132 {
133 assert((size_t)i <= columnEntries.size());
134 map<int,double>::iterator pos = columnEntries[i].find(j);
135 if (pos != columnEntries[i].end())
136 pos->second += value;
137 else
138 {
139 pair<int,double> entry(j,value);
140 columnEntries[i].insert(entry);
141 }
142 }
143
144 // remove the entry in row i, column j
RemoveEntry(int i,int j)145 void SparseMatrixOutline::RemoveEntry(int i, int j)
146 {
147 assert((size_t)i <= columnEntries.size());
148 columnEntries[i].erase(j);
149 }
150
151 // add a block matrix, starting at row i, and column j
AddBlockMatrix(int iStart,int jStart,const SparseMatrix * block,double scalarFactor)152 void SparseMatrixOutline::AddBlockMatrix(int iStart, int jStart, const SparseMatrix * block, double scalarFactor)
153 {
154 int nBlock = block->GetNumRows();
155 for(int i=0; i<nBlock; i++)
156 {
157 int rowLength = block->GetRowLength(i);
158 for(int j=0; j<rowLength; j++)
159 AddEntry(iStart + i, jStart + block->GetColumnIndex(i,j), scalarFactor * block->GetEntry(i,j));
160 }
161 }
162
MultiplyRow(int row,double scalar)163 void SparseMatrixOutline::MultiplyRow(int row, double scalar)
164 {
165 for(map<int,double>::iterator iter = columnEntries[row].begin(); iter != columnEntries[row].end(); iter++)
166 iter->second *= scalar;
167 }
168
AddBlock3x3Entry(int i,int j,const double * matrix3x3)169 void SparseMatrixOutline::AddBlock3x3Entry(int i, int j, const double * matrix3x3)
170 {
171 for(int k=0; k<3; k++)
172 for(int l=0; l<3; l++)
173 AddEntry(3*i+k,3*j+l,matrix3x3[3*k+l]);
174 }
175
AddBlock3x3Entry(int i,int j,double value)176 void SparseMatrixOutline::AddBlock3x3Entry(int i, int j, double value)
177 {
178 for(int k=0; k<3; k++)
179 for(int l=0; l<3; l++)
180 AddEntry(3*i+k,3*j+l,value);
181 }
182
GetEntry(int i,int j) const183 double SparseMatrixOutline::GetEntry(int i, int j) const
184 {
185 map<int,double>::const_iterator pos = columnEntries[i].find(j);
186 if (pos != columnEntries[i].end())
187 return (pos->second);
188 else
189 return 0;
190 }
191
GetNumColumns() const192 int SparseMatrixOutline::GetNumColumns() const
193 {
194 int numColumns = -1;
195 for(int i=0; i<numRows; i++)
196 {
197 map<int,double>::const_iterator j1;
198 // traverse all row entries
199 for(j1=columnEntries[i].begin(); j1 != columnEntries[i].end(); j1++)
200 if (j1->first > numColumns)
201 numColumns = j1->first;
202 }
203 return numColumns + 1;
204 }
205
Save(const char * filename,int oneIndexed) const206 int SparseMatrixOutline::Save(const char * filename, int oneIndexed) const
207 {
208 FILE * fout = fopen(filename, "w");
209 if (!fout)
210 return 1;
211
212 fprintf(fout, "%d\n%d\n", numRows, GetNumColumns());
213 for(int i=0; i<numRows; i++)
214 {
215 map<int,double>::const_iterator j1;
216 // traverse all row entries
217 for(j1=columnEntries[i].begin(); j1 != columnEntries[i].end(); ++j1)
218 fprintf(fout,"%d %d %.15f\n",i,j1->first + oneIndexed, j1->second + oneIndexed);
219 }
220 fclose(fout);
221
222 return 0;
223 }
224
Print() const225 void SparseMatrixOutline::Print() const
226 {
227 for (int i=0; i<numRows; i++)
228 {
229 for (int j=0; j<numRows; j++)
230 printf("%f ",GetEntry(i,j));
231 printf("\n");
232 }
233 }
234
GetNumEntries() const235 int SparseMatrixOutline::GetNumEntries() const
236 {
237 int num = 0;
238 for(int i=0; i<numRows; i++)
239 num += (int)(columnEntries[i].size());
240 return num;
241 }
242
SparseMatrix(const char * filename)243 SparseMatrix::SparseMatrix(const char * filename)
244 {
245 SparseMatrixOutline sparseMatrixOutline(filename);
246 InitFromOutline(&sparseMatrixOutline);
247 }
248
SparseMatrix(SparseMatrixOutline * sparseMatrixOutline)249 SparseMatrix::SparseMatrix(SparseMatrixOutline * sparseMatrixOutline)
250 {
251 InitFromOutline(sparseMatrixOutline);
252 }
253
254 // construct matrix from the outline
InitFromOutline(SparseMatrixOutline * sparseMatrixOutline)255 void SparseMatrix::InitFromOutline(SparseMatrixOutline * sparseMatrixOutline)
256 {
257 numRows = sparseMatrixOutline->GetNumRows();
258 Allocate();
259
260 for(int i=0; i<numRows; i++)
261 {
262 rowLength[i] = (int)(sparseMatrixOutline->columnEntries[i].size());
263 columnIndices[i] = (int*) malloc (sizeof(int) * rowLength[i]);
264 columnEntries[i] = (double*) malloc (sizeof(double) * rowLength[i]);
265
266 map<int,double>::iterator pos;
267 int j = 0;
268 int prev = -1;
269 for(pos = sparseMatrixOutline->columnEntries[i].begin(); pos != sparseMatrixOutline->columnEntries[i].end(); pos++)
270 {
271 columnIndices[i][j] = pos->first;
272 if (columnIndices[i][j] <= prev)
273 printf("Warning: entries not sorted in a row in a sparse matrix.\n");
274 prev = columnIndices[i][j];
275 columnEntries[i][j] = pos->second;
276 j++;
277 }
278 }
279 }
280
281 // construct matrix by specifying all entries: number of rows, length of each row, indices of columns of non-zero entries in each row, values of non-zero entries in each row
SparseMatrix(int numRows_,int * rowLength_,int ** columnIndices_,double ** columnEntries_,int shallowCopy)282 SparseMatrix::SparseMatrix(int numRows_, int * rowLength_, int ** columnIndices_, double ** columnEntries_, int shallowCopy) : numRows(numRows_)
283 {
284 if (shallowCopy)
285 {
286 rowLength = rowLength_;
287 columnIndices = columnIndices_;
288 columnEntries = columnEntries_;
289 numSubMatrixIDs = 0;
290 subMatrixIndices = NULL;
291 subMatrixIndexLengths = NULL;
292 subMatrixStartRow = NULL;
293 subMatrixNumRows = NULL;
294 superMatrixIndices = NULL;
295 superRows = NULL;
296 diagonalIndices = NULL;
297 transposedIndices = NULL;
298 return;
299 }
300
301 Allocate();
302
303 for (int i = 0; i<numRows; i++)
304 {
305 rowLength[i] = rowLength_[i];
306 columnIndices[i] = (int*) malloc(sizeof(int) * rowLength[i]);
307 columnEntries[i] = (double*) malloc(sizeof(double) * rowLength[i]);
308 for (int j = 0; j < rowLength[i]; j++)
309 {
310 columnIndices[i][j] = columnIndices_[i][j];
311 columnEntries[i][j] = columnEntries_[i][j];
312 }
313 }
314 }
315
316 // allocator
Allocate()317 void SparseMatrix::Allocate()
318 {
319 // compressed row storage
320 rowLength = (int*) malloc(sizeof(int) * numRows);
321 columnIndices = (int**) malloc(sizeof(int*) * numRows);
322 columnEntries = (double**) malloc(sizeof(double*) * numRows);
323 numSubMatrixIDs = 0;
324 subMatrixIndices = NULL;
325 subMatrixIndexLengths = NULL;
326 subMatrixStartRow = NULL;
327 subMatrixNumRows = NULL;
328 superMatrixIndices = NULL;
329 superRows = NULL;
330 diagonalIndices = NULL;
331 transposedIndices = NULL;
332 }
333
334 // destructor
~SparseMatrix()335 SparseMatrix::~SparseMatrix()
336 {
337 for(int i=0; i<numRows; i++)
338 {
339 free(columnIndices[i]);
340 free(columnEntries[i]);
341 }
342
343 if (subMatrixIndices != NULL)
344 {
345 for(int i=numSubMatrixIDs-1; i>=0; i--)
346 FreeSubMatrixIndices(i);
347 }
348
349 if (superRows != NULL)
350 {
351 for(int i=0; i<numRows; i++)
352 free(superMatrixIndices[i]);
353 free(superMatrixIndices);
354 free(superRows);
355 }
356
357 free(rowLength);
358 free(columnIndices);
359 free(columnEntries);
360 free(diagonalIndices);
361 FreeTranspositionIndices();
362 }
363
364 // copy constructor
SparseMatrix(const SparseMatrix & source)365 SparseMatrix::SparseMatrix(const SparseMatrix & source)
366 {
367 //printf("Copy constructor\n");fflush(NULL);
368 numRows = source.GetNumRows();
369
370 // compressed row storage
371 rowLength = (int*) malloc(sizeof(int) * numRows);
372 columnIndices = (int**) malloc(sizeof(int*) * numRows);
373 columnEntries = (double**) malloc(sizeof(double*) * numRows);
374
375 for(int i=0; i<numRows; i++)
376 {
377 rowLength[i] = source.rowLength[i];
378 columnIndices[i] = (int*) malloc (sizeof(int) * rowLength[i]);
379 columnEntries[i] = (double*) malloc (sizeof(double) * rowLength[i]);
380
381 for(int j=0; j < rowLength[i]; j++)
382 {
383 columnIndices[i][j] = source.columnIndices[i][j];
384 columnEntries[i][j] = source.columnEntries[i][j];
385 }
386 }
387
388 subMatrixIndices = NULL;
389 subMatrixIndexLengths = NULL;
390 subMatrixStartRow = NULL;
391 subMatrixNumRows = NULL;
392 numSubMatrixIDs = source.numSubMatrixIDs;
393 if (source.subMatrixIndices != NULL)
394 {
395 subMatrixIndices = (int***) malloc (sizeof(int**) * numSubMatrixIDs);
396 subMatrixIndexLengths = (int**) malloc (sizeof(int*) * numSubMatrixIDs);
397 subMatrixStartRow = (int*) malloc(sizeof(int) * numSubMatrixIDs);
398 subMatrixNumRows = (int*) malloc(sizeof(int) * numSubMatrixIDs);
399 memcpy(subMatrixStartRow, source.subMatrixStartRow, sizeof(int) * numSubMatrixIDs);
400 memcpy(subMatrixNumRows, source.subMatrixNumRows, sizeof(int) * numSubMatrixIDs);
401
402 for(int matrixID=0; matrixID < numSubMatrixIDs; matrixID++)
403 {
404 if (source.subMatrixIndices[matrixID] == NULL)
405 {
406 subMatrixIndices[matrixID] = NULL;
407 subMatrixIndexLengths[matrixID] = NULL;
408 continue;
409 }
410
411 subMatrixIndices[matrixID] = (int**) malloc(sizeof(int*) * numRows);
412 subMatrixIndexLengths[matrixID] = (int*) malloc(sizeof(int) * numRows);
413
414 for(int i=0; i<subMatrixNumRows[matrixID]; i++)
415 {
416 subMatrixIndexLengths[matrixID][i] = source.subMatrixIndexLengths[matrixID][i];
417 subMatrixIndices[matrixID][i] = (int*) malloc(sizeof(int) * subMatrixIndexLengths[matrixID][i]);
418 for(int j=0; j < subMatrixIndexLengths[matrixID][i]; j++)
419 subMatrixIndices[matrixID][i][j] = source.subMatrixIndices[matrixID][i][j];
420 }
421 }
422 }
423
424 superRows = NULL;
425 superMatrixIndices = NULL;
426 if (source.superRows != NULL)
427 {
428 superRows = (int*) malloc(sizeof(int) * numRows);
429 superMatrixIndices = (int**) malloc(sizeof(int*) * numRows);
430 for(int i=0; i<numRows; i++)
431 {
432 superRows[i] = source.superRows[i];
433 superMatrixIndices[i] = (int*) malloc(sizeof(int) * rowLength[i]);
434 for(int j=0; j < rowLength[i]; j++)
435 {
436 superMatrixIndices[i][j] = source.superMatrixIndices[i][j];
437 }
438 }
439 }
440
441 diagonalIndices = NULL;
442 if (source.diagonalIndices != NULL)
443 {
444 diagonalIndices = (int*) malloc (sizeof(int) * numRows);
445 memcpy(diagonalIndices, source.diagonalIndices, sizeof(int) * numRows);
446 }
447
448 transposedIndices = NULL;
449 if (source.transposedIndices != NULL)
450 {
451 transposedIndices = (int**) malloc (sizeof(int*) * numRows);
452 for(int i=0; i<numRows; i++)
453 {
454 transposedIndices[i] = (int*) malloc (sizeof(int) * rowLength[i]);
455 for(int j=0; j<rowLength[i]; j++)
456 transposedIndices[i][j] = source.transposedIndices[i][j];
457 }
458 }
459 }
460
MultiplyVector(int startRow,int endRow,const double * vector,double * result) const461 void SparseMatrix::MultiplyVector(int startRow, int endRow, const double * vector, double * result) const // result = A(startRow:endRow-1,:) * vector
462 {
463 for(int i=startRow; i<endRow; i++)
464 {
465 result[i-startRow] = 0;
466 for(int j=0; j < rowLength[i]; j++)
467 result[i-startRow] += vector[columnIndices[i][j]] * columnEntries[i][j];
468 }
469 }
470
MultiplyVector(const double * vector,double * result) const471 void SparseMatrix::MultiplyVector(const double * vector, double * result) const
472 {
473 for(int i=0; i<numRows; i++)
474 {
475 result[i] = 0;
476 for(int j=0; j < rowLength[i]; j++)
477 result[i] += vector[columnIndices[i][j]] * columnEntries[i][j];
478 }
479 }
480
MultiplyVectorAdd(const double * vector,double * result) const481 void SparseMatrix::MultiplyVectorAdd(const double * vector, double * result) const
482 {
483 for(int i=0; i<numRows; i++)
484 for(int j=0; j < rowLength[i]; j++)
485 result[i] += vector[columnIndices[i][j]] * columnEntries[i][j];
486 }
487
TransposeMultiplyVector(const double * vector,int resultLength,double * result) const488 void SparseMatrix::TransposeMultiplyVector(const double * vector, int resultLength, double * result) const
489 {
490 for(int i=0; i<resultLength; i++)
491 result[i] = 0;
492
493 for(int i=0; i<numRows; i++)
494 for(int j=0; j < rowLength[i]; j++)
495 result[columnIndices[i][j]] += vector[i] * columnEntries[i][j];
496 }
497
TransposeMultiplyVectorAdd(const double * vector,double * result) const498 void SparseMatrix::TransposeMultiplyVectorAdd(const double * vector, double * result) const
499 {
500 for(int i=0; i<numRows; i++)
501 {
502 for(int j=0; j < rowLength[i]; j++)
503 {
504 result[columnIndices[i][j]] += vector[i] * columnEntries[i][j];
505 }
506 }
507 }
508
MultiplyMatrix(int numDenseRows,int numDenseColumns,const double * denseMatrix,double * result) const509 void SparseMatrix::MultiplyMatrix(int numDenseRows, int numDenseColumns, const double * denseMatrix, double * result) const
510 {
511 for(int column=0; column<numDenseColumns; column++)
512 MultiplyVector(&denseMatrix[numDenseRows * column], &result[numRows * column]);
513 }
514
MultiplyMatrixAdd(int numDenseRows,int numDenseColumns,const double * denseMatrix,double * result) const515 void SparseMatrix::MultiplyMatrixAdd(int numDenseRows, int numDenseColumns, const double * denseMatrix, double * result) const
516 {
517 for(int column=0; column<numDenseColumns; column++)
518 MultiplyVectorAdd(&denseMatrix[numDenseRows * column], &result[numRows * column]);
519 }
520
521 // result = A * trans(denseMatrix)
522 // trans(denseMatrix) is a dense matrix with 'numDenseColumns' columns, result is a numRows x numDenseColumns dense matrix
MultiplyMatrixTranspose(int numDenseColumns,const double * denseMatrix,double * result) const523 void SparseMatrix::MultiplyMatrixTranspose(int numDenseColumns, const double * denseMatrix, double * result) const
524 {
525 memset(result, 0, sizeof(double) * numRows * numDenseColumns);
526 for(int column=0; column<numDenseColumns; column++)
527 for(int i=0; i<numRows; i++)
528 for(int j=0; j < rowLength[i]; j++)
529 result[numRows * column + i] += denseMatrix[numDenseColumns * columnIndices[i][j] + column] * columnEntries[i][j];
530 }
531
QuadraticForm(const double * vector) const532 double SparseMatrix::QuadraticForm(const double * vector) const
533 {
534 double result = 0;
535 for(int i=0; i<numRows; i++)
536 {
537 for(int j=0; j < rowLength[i]; j++)
538 {
539 int index = columnIndices[i][j];
540 if (index < i)
541 continue;
542 if (index == i)
543 result += columnEntries[i][j] * vector[i] * vector[index];
544 else
545 result += 2.0 * columnEntries[i][j] * vector[i] * vector[index];
546 }
547 }
548
549 return result;
550 }
551
NormalizeVector(double * vector) const552 void SparseMatrix::NormalizeVector(double * vector) const
553 {
554 double norm = sqrt(QuadraticForm(vector));
555 for(int i=0; i<numRows; i++)
556 vector[i] /= norm;
557 }
558
operator +(const SparseMatrix & mat2) const559 SparseMatrix SparseMatrix::operator+(const SparseMatrix & mat2) const
560 {
561 SparseMatrix result(*this);
562 for(int i=0; i<numRows; i++)
563 for(int j=0; j < rowLength[i]; j++)
564 result.columnEntries[i][j] += mat2.columnEntries[i][j];
565 return result;
566 }
567
operator -(const SparseMatrix & mat2) const568 SparseMatrix SparseMatrix::operator-(const SparseMatrix & mat2) const
569 {
570 SparseMatrix result(*this);
571 for(int i=0; i<numRows; i++)
572 for(int j=0; j < rowLength[i]; j++)
573 result.columnEntries[i][j] -= mat2.columnEntries[i][j];
574 return result;
575 }
576
operator *(const double alpha,const SparseMatrix & mat2)577 SparseMatrix operator* (const double alpha, const SparseMatrix & mat2)
578 {
579 SparseMatrix result(mat2);
580 for(int i=0; i<result.numRows; i++)
581 for(int j=0; j < result.rowLength[i]; j++)
582 result.columnEntries[i][j] *= alpha;
583 return result;
584 }
585
operator *=(const double alpha)586 SparseMatrix & SparseMatrix::operator*=(const double alpha)
587 {
588 for(int i=0; i<numRows; i++)
589 for(int j=0; j < rowLength[i]; j++)
590 columnEntries[i][j] *= alpha;
591 return *this;
592 }
593
operator +=(const SparseMatrix & mat2)594 SparseMatrix & SparseMatrix::operator+=(const SparseMatrix & mat2)
595 {
596 for(int i=0; i<numRows; i++)
597 for(int j=0; j < rowLength[i]; j++)
598 columnEntries[i][j] += mat2.columnEntries[i][j];
599 return *this;
600 }
601
operator -=(const SparseMatrix & mat2)602 SparseMatrix & SparseMatrix::operator-=(const SparseMatrix & mat2)
603 {
604 for(int i=0; i<numRows; i++)
605 for(int j=0; j < rowLength[i]; j++)
606 columnEntries[i][j] -= mat2.columnEntries[i][j];
607 return *this;
608 }
609
operator =(const SparseMatrix & source)610 SparseMatrix & SparseMatrix::operator=(const SparseMatrix & source)
611 {
612 for(int i=0; i<numRows; i++)
613 {
614 for(int j=0; j < rowLength[i]; j++)
615 columnEntries[i][j] = source.columnEntries[i][j];
616 }
617
618 return *this;
619 }
620
operator ==(const SparseMatrix & mat2) const621 bool SparseMatrix::operator==(const SparseMatrix & mat2) const
622 {
623 if(numRows != mat2.numRows)
624 return false;
625 for(int i = 0; i < numRows; i++)
626 {
627 if(rowLength[i] != mat2.rowLength[i])
628 return false;
629 for(int j = 0; j < rowLength[i]; j++)
630 {
631 if(columnIndices[i][j] != mat2.columnIndices[i][j])
632 return false;
633 if(columnEntries[i][j] != mat2.columnEntries[i][j])
634 return false;
635 }
636 }
637 return true;
638 }
639
SameStructure(const SparseMatrix & mat2) const640 bool SparseMatrix::SameStructure(const SparseMatrix & mat2) const
641 {
642 if(numRows != mat2.numRows)
643 return false;
644 for(int i = 0; i < numRows; i++)
645 {
646 if(rowLength[i] != mat2.rowLength[i])
647 return false;
648 for(int j = 0; j < rowLength[i]; j++)
649 {
650 if(columnIndices[i][j] != mat2.columnIndices[i][j])
651 return false;
652 }
653 }
654 return true;
655 }
656
ScalarMultiply(const double alpha,SparseMatrix * dest)657 void SparseMatrix::ScalarMultiply(const double alpha, SparseMatrix * dest)
658 {
659 if (dest == NULL)
660 dest = this;
661
662 for(int i=0; i<numRows; i++)
663 for(int j=0; j < rowLength[i]; j++)
664 dest->columnEntries[i][j] = columnEntries[i][j] * alpha;
665 }
666
ScalarMultiplyAdd(const double alpha,SparseMatrix * dest)667 void SparseMatrix::ScalarMultiplyAdd(const double alpha, SparseMatrix * dest)
668 {
669 if (dest == NULL)
670 dest = this;
671
672 for(int i=0; i<numRows; i++)
673 for(int j=0; j < rowLength[i]; j++)
674 dest->columnEntries[i][j] += columnEntries[i][j] * alpha;
675 }
676
ResetToZero()677 void SparseMatrix::ResetToZero()
678 {
679 for(int i=0; i<numRows; i++)
680 memset(columnEntries[i], 0, sizeof(double) * rowLength[i]);
681 }
682
ResetRowToZero(int row)683 void SparseMatrix::ResetRowToZero(int row)
684 {
685 memset(columnEntries[row], 0, sizeof(double) * rowLength[row]);
686 }
687
Print(int sparsePrint) const688 void SparseMatrix::Print(int sparsePrint) const
689 {
690 if (sparsePrint)
691 {
692 for (int i=0; i<numRows; i++)
693 for(int j=0; j< rowLength[i]; j++)
694 printf("%d %d %G\n", i, columnIndices[i][j], columnEntries[i][j]);
695 }
696 else
697 {
698 int numColumns = GetNumColumns();
699 for (int i=0; i<numRows; i++)
700 {
701 int index = 0;
702 for(int j=0; j< rowLength[i]; j++)
703 {
704 while (index < columnIndices[i][j])
705 {
706 index++;
707 printf("%f,",0.0);
708 }
709 printf("%f,",columnEntries[i][j]);
710 index++;
711 }
712
713 while (index < numColumns)
714 {
715 index++;
716 printf("%f,",0.0);
717 }
718
719 printf("\n");
720 }
721 }
722 }
723
724 // seeks the element in column jDense (in row "row")
725 // if not found, returns -1
GetInverseIndex(int row,int jDense) const726 int SparseMatrix::GetInverseIndex(int row, int jDense) const
727 {
728 for(int j=0; j < rowLength[row]; j++)
729 if (columnIndices[row][j] == jDense)
730 return j;
731
732 return -1;
733 }
734
BuildDiagonalIndices()735 void SparseMatrix::BuildDiagonalIndices()
736 {
737 if (diagonalIndices != NULL)
738 return;
739
740 diagonalIndices = (int*) malloc (sizeof(int) * numRows);
741 for(int i=0; i<numRows; i++)
742 diagonalIndices[i] = GetInverseIndex(i,i);
743 }
744
FreeDiagonalIndices()745 void SparseMatrix::FreeDiagonalIndices()
746 {
747 free(diagonalIndices);
748 }
749
GetDiagonal(double * diagonal) const750 void SparseMatrix::GetDiagonal(double * diagonal) const
751 {
752 if (diagonalIndices != NULL)
753 {
754 for(int i=0; i<numRows; i++)
755 diagonal[i] = columnEntries[i][diagonalIndices[i]];
756 }
757 else
758 {
759 for(int i=0; i<numRows; i++)
760 for(int j=0; j<GetRowLength(i); j++)
761 {
762 if (GetColumnIndex(i, j) == i)
763 diagonal[i] = columnEntries[i][j];
764 }
765 }
766 }
767
AddDiagonalMatrix(double * diagonalMatrix)768 void SparseMatrix::AddDiagonalMatrix(double * diagonalMatrix)
769 {
770 if (diagonalIndices != NULL)
771 {
772 for(int i=0; i<numRows; i++)
773 columnEntries[i][diagonalIndices[i]] += diagonalMatrix[i];
774 }
775 else
776 {
777 for(int i=0; i<numRows; i++)
778 for(int j=0; j<GetRowLength(i); j++)
779 {
780 if (GetColumnIndex(i, j) == i)
781 columnEntries[i][j] += diagonalMatrix[i];
782 }
783 }
784 }
785
AddDiagonalMatrix(double constDiagonalElement)786 void SparseMatrix::AddDiagonalMatrix(double constDiagonalElement)
787 {
788 if (diagonalIndices != NULL)
789 {
790 for(int i=0; i<numRows; i++)
791 columnEntries[i][diagonalIndices[i]] += constDiagonalElement;
792 }
793 else
794 {
795 for(int i=0; i<numRows; i++)
796 for(int j=0; j<GetRowLength(i); j++)
797 {
798 if (GetColumnIndex(i, j) == i)
799 columnEntries[i][j] += constDiagonalElement;
800 }
801 }
802 }
803
GetNumEntries() const804 int SparseMatrix::GetNumEntries() const
805 {
806 int num = 0;
807 for(int i=0; i<numRows; i++)
808 num += rowLength[i];
809
810 return num;
811 }
812
SumEntries() const813 double SparseMatrix::SumEntries() const
814 {
815 double sum=0;
816 for(int i=0; i<numRows; i++)
817 for(int j=0; j<rowLength[i]; j++)
818 sum += columnEntries[i][j];
819
820 return sum;
821 }
822
SumRowEntries(double * rowSums) const823 void SparseMatrix::SumRowEntries(double * rowSums) const
824 {
825 for(int i=0; i<numRows; i++)
826 {
827 double sum=0;
828 for(int j=0; j<rowLength[i]; j++)
829 sum += columnEntries[i][j];
830 rowSums[i] = sum;
831 }
832 }
833
MakeLinearDataArray(double * data) const834 void SparseMatrix::MakeLinearDataArray(double * data) const
835 {
836 int count=0;
837 for(int i=0; i<numRows; i++)
838 {
839 for(int j=0; j<rowLength[i]; j++)
840 {
841 data[count] = columnEntries[i][j];
842 count++;
843 }
844 }
845 }
846
MakeLinearRowIndexArray(double * indices) const847 void SparseMatrix::MakeLinearRowIndexArray(double * indices) const
848 {
849 int count=0;
850 for(int i=0; i<numRows; i++)
851 {
852 for(int j=0; j<rowLength[i]; j++)
853 {
854 indices[count] = i;
855 count++;
856 }
857 }
858 }
859
MakeLinearColumnIndexArray(double * indices) const860 void SparseMatrix::MakeLinearColumnIndexArray(double * indices) const
861 {
862 int count=0;
863 for(int i=0; i<numRows; i++)
864 {
865 for(int j=0; j<rowLength[i]; j++)
866 {
867 indices[count] = columnIndices[i][j];
868 count++;
869 }
870 }
871 }
872
MakeLinearRowIndexArray(int * indices) const873 void SparseMatrix::MakeLinearRowIndexArray(int * indices) const
874 {
875 int count=0;
876 for(int i=0; i<numRows; i++)
877 {
878 for(int j=0; j<rowLength[i]; j++)
879 {
880 indices[count] = i;
881 count++;
882 }
883 }
884 }
885
MakeLinearColumnIndexArray(int * indices) const886 void SparseMatrix::MakeLinearColumnIndexArray(int * indices) const
887 {
888 int count=0;
889 for(int i=0; i<numRows; i++)
890 {
891 for(int j=0; j<rowLength[i]; j++)
892 {
893 indices[count] = columnIndices[i][j];
894 count++;
895 }
896 }
897 }
898
FreeTranspositionIndices()899 void SparseMatrix::FreeTranspositionIndices()
900 {
901 if (transposedIndices == NULL)
902 return;
903
904 for(int i=0; i<numRows; i++)
905 free(transposedIndices[i]);
906 free(transposedIndices);
907
908 transposedIndices = NULL;
909 }
910
BuildTranspositionIndices()911 void SparseMatrix::BuildTranspositionIndices()
912 {
913 if (transposedIndices != NULL)
914 return;
915
916 transposedIndices = (int**) malloc (sizeof(int*) * numRows);
917 int * buffer = (int*) calloc (GetNumColumns(), sizeof(int));
918
919 for(int i=0; i<numRows; i++)
920 {
921 transposedIndices[i] = (int*) malloc (sizeof(int) * rowLength[i]);
922 for(int j=0; j<rowLength[i]; j++)
923 {
924 transposedIndices[i][j] = buffer[columnIndices[i][j]];
925 buffer[columnIndices[i][j]]++;
926 }
927 }
928
929 free(buffer);
930 }
931
SkewSymmetricCheck()932 double SparseMatrix::SkewSymmetricCheck()
933 {
934 double maxEntry = 0;
935
936 BuildTranspositionIndices();
937
938 for(int i=0; i<numRows; i++)
939 {
940 for(int j=0; j<GetRowLength(i); j++)
941 {
942 double entry1 = GetEntry(i, j);
943 int tindex = TransposedIndex(i, j);
944 double entry2 = GetEntry(GetColumnIndex(i,j), tindex);
945
946 // entry1 + entry2 should be zero
947 if (fabs(entry1 + entry2) > maxEntry)
948 maxEntry = fabs(entry1 + entry2);
949 }
950 }
951
952 return maxEntry;
953 }
954
SymmetrizeMatrix()955 void SparseMatrix::SymmetrizeMatrix()
956 {
957 BuildTranspositionIndices();
958
959 for(int i=0; i<numRows; i++)
960 {
961 for(int j=0; j<rowLength[i]; j++)
962 {
963 int jAbs = columnIndices[i][j];
964
965 if (jAbs >= i)
966 break;
967
968 // copy elt (jAbs,i) into position (i,jAbs)
969 columnEntries[i][j] = columnEntries[jAbs][TransposedIndex(i,j)];
970 }
971 }
972 }
973
GetMaxAbsEntry() const974 double SparseMatrix::GetMaxAbsEntry() const
975 {
976 double max = 0;
977 for(int i=0; i<numRows; i++)
978 {
979 for(int j=0; j<rowLength[i]; j++)
980 {
981 double el = fabs(GetEntry(i,j));
982 if (el > max)
983 max = el;
984 }
985 }
986
987 return max;
988 }
989
GetRowNorm2(int row) const990 double SparseMatrix::GetRowNorm2(int row) const
991 {
992 double norm2 = 0;
993 for(int j=0; j<rowLength[row]; j++)
994 {
995 double el = columnEntries[row][j];
996 norm2 += el * el;
997 }
998 return norm2;
999 }
1000
1001 // solve M * x = b
1002 // ASSUMES the sparse matrix is diagonal !!!
DiagonalSolve(double * rhs) const1003 void SparseMatrix::DiagonalSolve(double * rhs) const
1004 {
1005 for(int i=0; i<numRows; i++)
1006 rhs[i] /= columnEntries[i][0]; // the diagonal element
1007 }
1008
BuildRenumberingVector(int nConstrained,int nSuper,int numFixedDOFs,int * fixedDOFs,int ** superDOFs,int oneIndexed)1009 void SparseMatrix::BuildRenumberingVector(int nConstrained, int nSuper, int numFixedDOFs, int * fixedDOFs, int ** superDOFs, int oneIndexed)
1010 {
1011 // superRows[i] is the row index in the super matrix corresponsing to row i of constrained matrix
1012 (*superDOFs) = (int*) malloc (sizeof(int) * nConstrained);
1013 int constrainedDOF = 0;
1014 int superDOF = 0;
1015 for(int i=0; i<numFixedDOFs; i++)
1016 {
1017 int nextSuperDOF = fixedDOFs[i];
1018 nextSuperDOF -= oneIndexed;
1019 if ( (nextSuperDOF >= nSuper) || (nextSuperDOF < 0) )
1020 {
1021 printf("Error: invalid fixed super DOF %d specified.\n", nextSuperDOF);
1022 exit(1);
1023 }
1024
1025 while (superDOF < nextSuperDOF)
1026 {
1027 if (constrainedDOF >= nConstrained)
1028 {
1029 printf("Error: too many DOFs specified.\n");
1030 exit(1);
1031 }
1032 (*superDOFs)[constrainedDOF] = superDOF;
1033 constrainedDOF++;
1034 superDOF++;
1035 }
1036
1037 superDOF++; // skip the deselected DOF
1038 }
1039 while (superDOF < nSuper)
1040 {
1041 if (constrainedDOF >= nConstrained)
1042 {
1043 printf("Error: too many DOFs specified.\n");
1044 exit(1);
1045 }
1046 (*superDOFs)[constrainedDOF] = superDOF;
1047
1048 constrainedDOF++;
1049 superDOF++;
1050 }
1051 }
1052
BuildSuperMatrixIndices(int numFixedRowColumns,int * fixedRowColumns,const SparseMatrix * superMatrix,int oneIndexed)1053 void SparseMatrix::BuildSuperMatrixIndices(int numFixedRowColumns, int * fixedRowColumns, const SparseMatrix * superMatrix, int oneIndexed)
1054 {
1055 BuildSuperMatrixIndices(numFixedRowColumns, fixedRowColumns, numFixedRowColumns, fixedRowColumns, superMatrix, oneIndexed);
1056 }
1057
BuildSuperMatrixIndices(int numFixedRows,int * fixedRows,int numFixedColumns,int * fixedColumns,const SparseMatrix * superMatrix,int oneIndexed)1058 void SparseMatrix::BuildSuperMatrixIndices(int numFixedRows, int * fixedRows, int numFixedColumns, int * fixedColumns, const SparseMatrix * superMatrix, int oneIndexed)
1059 {
1060 int numSuperColumns = superMatrix->GetNumColumns();
1061 int numColumns = numSuperColumns - numFixedColumns;
1062
1063 if ((numRows + numFixedRows != superMatrix->numRows) || (GetNumColumns() + numFixedColumns > numSuperColumns) )
1064 {
1065 printf("Error in BuildSuperMatrixIndices: number of constrained DOFs does not match the size of the two matrices.\n");
1066 printf("num rows: %d num fixed rows in super matrix: %d num rows in super matrix: %d\n", numRows, numFixedRows, superMatrix->numRows);
1067 printf("num columns: %d num fixed columns in super matrix: %d num columns in super matrix: %d\n", numColumns, numFixedColumns, numSuperColumns);
1068 exit(1);
1069 }
1070
1071 // build row renumbering function:
1072 BuildRenumberingVector(numRows, superMatrix->numRows, numFixedRows, fixedRows, &superRows, oneIndexed);
1073 // build column renumbering function:
1074 int * superColumns_;
1075 BuildRenumberingVector(numColumns, numSuperColumns, numFixedColumns, fixedColumns, &superColumns_, oneIndexed);
1076
1077 // superRows[i] is the row index in the super matrix corresponsing to row i of constrained matrix
1078 // superColumns_[i] is the dense column index in the super matrix corresponsing to the dense column i of constrained matrix
1079
1080 // build column indices
1081 superMatrixIndices = (int**) malloc (sizeof(int*) * numRows);
1082 for(int i=0; i < numRows; i++)
1083 {
1084 superMatrixIndices[i] = (int*) malloc (sizeof(int) * rowLength[i]);
1085 for(int j=0; j < rowLength[i]; j++)
1086 {
1087 int iConstrained = i;
1088 int jConstrainedDense = columnIndices[iConstrained][j];
1089 int iSuper = superRows[iConstrained];
1090 int jSuperDense = superColumns_[jConstrainedDense];
1091 int jSuper = superMatrix->GetInverseIndex(iSuper, jSuperDense);
1092 if (jSuper < 0)
1093 {
1094 printf("Error in BuildSuperMatrixIndices: failed to compute inverse index.\n");
1095 printf("i=%d j=%d iConstrained=%d jConstrainedDense=%d iSuper=%d jSuperDense=%d jSuper=%d\n", i, j, iConstrained, jConstrainedDense, iSuper, jSuperDense, jSuper);
1096 fflush(NULL);
1097 exit(1);
1098 }
1099 superMatrixIndices[i][j] = jSuper;
1100 }
1101 }
1102
1103 free(superColumns_);
1104 }
1105
AssignSuperMatrix(const SparseMatrix & superMatrix)1106 void SparseMatrix::AssignSuperMatrix(const SparseMatrix & superMatrix)
1107 {
1108 for(int i=0; i<numRows; i++)
1109 {
1110 double * row = superMatrix.columnEntries[superRows[i]];
1111 int * indices = superMatrixIndices[i];
1112 for(int j=0; j < rowLength[i]; j++)
1113 columnEntries[i][j] = row[indices[j]];
1114 }
1115 }
1116
BuildSubMatrixIndices(const SparseMatrix & submatrix,int subMatrixID,int startRow,int startColumn)1117 void SparseMatrix::BuildSubMatrixIndices(const SparseMatrix & submatrix, int subMatrixID, int startRow, int startColumn)
1118 {
1119 if (subMatrixID >= numSubMatrixIDs)
1120 {
1121 subMatrixIndices = (int***) realloc (subMatrixIndices, sizeof(int**) * (subMatrixID + 1));
1122 subMatrixIndexLengths = (int**) realloc (subMatrixIndexLengths, sizeof(int*) * (subMatrixID + 1));
1123 subMatrixStartRow = (int*) realloc(subMatrixStartRow, sizeof(int) * (subMatrixID + 1));
1124 subMatrixNumRows = (int*) realloc(subMatrixNumRows, sizeof(int) * (subMatrixID + 1));
1125 for(int i=numSubMatrixIDs; i <= subMatrixID; i++)
1126 {
1127 subMatrixIndices[i] = NULL;
1128 subMatrixIndexLengths[i] = NULL;
1129 subMatrixStartRow[i] = 0;
1130 subMatrixNumRows[i] = 0;
1131 }
1132 numSubMatrixIDs = subMatrixID + 1;
1133 }
1134 else if ((subMatrixIndices[subMatrixID] != NULL) || (subMatrixIndexLengths[subMatrixID] != NULL))
1135 {
1136 free(subMatrixIndices[subMatrixID]);
1137 free(subMatrixIndexLengths[subMatrixID]);
1138 subMatrixIndices[subMatrixID] = NULL;
1139 subMatrixIndexLengths[subMatrixID] = NULL;
1140 subMatrixStartRow[subMatrixID] = 0;
1141 subMatrixNumRows[subMatrixID] = 0;
1142 //printf("Warning: old submatrix indices (matrixID %d) have not been de-allocated.\n", subMatrixID);
1143 }
1144
1145 subMatrixStartRow[subMatrixID] = startRow;
1146 subMatrixNumRows[subMatrixID] = submatrix.numRows;
1147 subMatrixIndices[subMatrixID] = (int**) malloc (sizeof(int*) * subMatrixNumRows[subMatrixID]);
1148 subMatrixIndexLengths[subMatrixID] = (int*) malloc (sizeof(int) * subMatrixNumRows[subMatrixID]);
1149
1150 int endRow = startRow + subMatrixNumRows[subMatrixID];
1151 if(endRow > numRows)
1152 {
1153 printf("Error (BuildSubMatrixIndices): given submatrix placed at startRow %d exceeds the length of this matrix: %d\n", startRow, numRows);
1154 exit(1);
1155 }
1156
1157 for(int i=0; i<subMatrixNumRows[subMatrixID]; i++)
1158 {
1159 //begin at the startRow, find correspondence to subMatrix at each row
1160 subMatrixIndices[subMatrixID][i] = (int*) malloc (sizeof(int) * submatrix.rowLength[i]);
1161 subMatrixIndexLengths[subMatrixID][i] = submatrix.rowLength[i];
1162 int * indices = submatrix.columnIndices[i];
1163
1164 // this assumes that the column indices are sorted (which they always are)
1165 int col = 0;
1166 int * colIndices = columnIndices[startRow + i];
1167
1168 for(int j=0; j < submatrix.rowLength[i]; j++)
1169 {
1170 int targetColumn = startColumn + indices[j];
1171 while (colIndices[col] < targetColumn)
1172 col++;
1173
1174 subMatrixIndices[subMatrixID][i][j] = (colIndices[col] == targetColumn) ? col : -1;
1175
1176 // finds the position in row i of element with column index jDense
1177 // int inverseIndex(int i, int jDense);
1178 //subMatrixIndices[subMatrixID][i][j] = GetInverseIndex(startRow + i, startColumn + indices[j]);
1179 if (subMatrixIndices[subMatrixID][i][j] == -1)
1180 {
1181 printf("Error (BuildSubMatrixIndices): given matrix is not a submatrix of this matrix. The following index does not exist in this matrix: (%d,%d)\n", startRow + i, startColumn + indices[j]);
1182 exit(1);
1183 }
1184 }
1185 }
1186 }
1187
FreeSubMatrixIndices(int subMatrixID)1188 void SparseMatrix::FreeSubMatrixIndices(int subMatrixID)
1189 {
1190 if (subMatrixID >= numSubMatrixIDs)
1191 {
1192 printf("Warning: attempted to free submatrix index that does not exist.\n");
1193 return;
1194 }
1195
1196 if (subMatrixIndices[subMatrixID] != NULL)
1197 {
1198 for(int i=0; i<subMatrixNumRows[subMatrixID]; i++)
1199 free(subMatrixIndices[subMatrixID][i]);
1200 free(subMatrixIndices[subMatrixID]);
1201 free(subMatrixIndexLengths[subMatrixID]);
1202 subMatrixIndices[subMatrixID] = NULL;
1203 subMatrixIndexLengths[subMatrixID] = NULL;
1204 subMatrixStartRow[subMatrixID] = 0;
1205 subMatrixNumRows[subMatrixID] = 0;
1206 }
1207
1208 // check if this was the largest index
1209 for(int i=numSubMatrixIDs-1; i>=0; i--)
1210 {
1211 if (subMatrixIndices[i] != NULL)
1212 {
1213 numSubMatrixIDs = i + 1;
1214 subMatrixIndices = (int***) realloc (subMatrixIndices, sizeof(int**) * numSubMatrixIDs);
1215 subMatrixIndexLengths = (int**) realloc (subMatrixIndexLengths, sizeof(int*) * numSubMatrixIDs);
1216 subMatrixStartRow = (int*) realloc (subMatrixStartRow, sizeof(int) * numSubMatrixIDs);
1217 subMatrixNumRows = (int*) realloc (subMatrixNumRows, sizeof(int) * numSubMatrixIDs);
1218 break;
1219 }
1220
1221 numSubMatrixIDs = i;
1222 }
1223
1224 if (numSubMatrixIDs == 0)
1225 {
1226 free(subMatrixIndices);
1227 free(subMatrixIndexLengths);
1228 free(subMatrixStartRow);
1229 free(subMatrixNumRows);
1230 subMatrixIndices = NULL;
1231 subMatrixIndexLengths = NULL;
1232 subMatrixStartRow = NULL;
1233 subMatrixNumRows = NULL;
1234 }
1235 }
1236
AssignSubMatrix(const SparseMatrix & submatrix,int subMatrixID)1237 void SparseMatrix::AssignSubMatrix(const SparseMatrix & submatrix, int subMatrixID)
1238 {
1239 int startRow = subMatrixStartRow[subMatrixID];
1240 for(int i=0; i<subMatrixNumRows[subMatrixID]; i++)
1241 {
1242 int * indices = subMatrixIndices[subMatrixID][i];
1243 for(int j=0; j < submatrix.rowLength[i]; j++)
1244 columnEntries[startRow + i][indices[j]] = submatrix.columnEntries[i][j];
1245 }
1246 }
1247
AssignToSubMatrix(SparseMatrix & submatrix,int subMatrixID) const1248 void SparseMatrix::AssignToSubMatrix(SparseMatrix & submatrix, int subMatrixID) const
1249 {
1250 int startRow = subMatrixStartRow[subMatrixID];
1251 for(int i=0; i<subMatrixNumRows[subMatrixID]; i++)
1252 {
1253 int * indices = subMatrixIndices[subMatrixID][i];
1254 for(int j=0; j < submatrix.rowLength[i]; j++)
1255 submatrix.columnEntries[i][j] = columnEntries[startRow + i][indices[j]];
1256 }
1257 }
1258
AddSubMatrix(double factor,SparseMatrix & submatrix,int subMatrixID)1259 SparseMatrix & SparseMatrix::AddSubMatrix(double factor, SparseMatrix & submatrix, int subMatrixID)
1260 {
1261 int startRow = subMatrixStartRow[subMatrixID];
1262 for(int i=0; i<subMatrixNumRows[subMatrixID]; i++)
1263 {
1264 int * indices = subMatrixIndices[subMatrixID][i];
1265 for(int j=0; j < submatrix.rowLength[i]; j++)
1266 columnEntries[startRow + i][indices[j]] += factor * submatrix.columnEntries[i][j];
1267 }
1268
1269 return *this;
1270 }
1271
GetNumLowerTriangleEntries() const1272 int SparseMatrix::GetNumLowerTriangleEntries() const
1273 {
1274 int num = 0;
1275 for(int i=0; i<numRows; i++)
1276 {
1277 for(int j=0; j < rowLength[i]; j++)
1278 {
1279 if (columnIndices[i][j] <= i)
1280 num++;
1281 }
1282 }
1283 return num;
1284 }
1285
GetNumUpperTriangleEntries() const1286 int SparseMatrix::GetNumUpperTriangleEntries() const
1287 {
1288 int num = 0;
1289 for(int i=0; i<numRows; i++)
1290 {
1291 for(int j=0; j < rowLength[i]; j++)
1292 {
1293 if (columnIndices[i][j] >= i)
1294 num++;
1295 }
1296 }
1297 return num;
1298 }
1299
GenerateNAGFormat(double * a,int * irow,int * icol,int * istr) const1300 int SparseMatrix::GenerateNAGFormat(double * a, int * irow, int * icol, int * istr) const
1301 {
1302 int num = 0;
1303 for(int i=0; i<numRows; i++)
1304 {
1305 istr[i] = num; // starting address of row i
1306 for(int j=0; j < rowLength[i]; j++)
1307 {
1308 if (columnIndices[i][j] <= i) // over lower triangle
1309 {
1310 a[num] = columnEntries[i][j];
1311 irow[num] = i+1; // NAG is 1-indexed
1312 icol[num] = columnIndices[i][j]+1; // NAG is 1-indexed
1313 num++;
1314 }
1315 }
1316 }
1317
1318 istr[numRows] = num;
1319
1320 return num;
1321 }
1322
GenerateCompressedRowMajorFormat(double * a,int * ia,int * ja,int upperTriangleOnly,int oneIndexed) const1323 void SparseMatrix::GenerateCompressedRowMajorFormat(double * a, int * ia, int * ja, int upperTriangleOnly, int oneIndexed) const
1324 {
1325 int count = 0;
1326 for(int row=0; row<numRows; row++)
1327 {
1328 if (ia != NULL)
1329 ia[row] = count + oneIndexed;
1330
1331 int rowLength = GetRowLength(row);
1332 for(int j=0; j< rowLength; j++)
1333 {
1334 if ((!upperTriangleOnly) || (columnIndices[row][j] >= row))
1335 {
1336 if (a != NULL)
1337 a[count] = columnEntries[row][j];
1338 if (ja != NULL)
1339 ja[count] = columnIndices[row][j] + oneIndexed;
1340 count++;
1341 }
1342 }
1343 }
1344
1345 if (ia != NULL)
1346 ia[numRows] = count + oneIndexed;
1347 }
1348
GenerateCompressedRowMajorFormat_four_array(double * values,int * columns,int * pointerB,int * pointerE,int upperTriangleOnly,int oneIndexed) const1349 void SparseMatrix::GenerateCompressedRowMajorFormat_four_array(double * values, int * columns, int * pointerB, int * pointerE, int upperTriangleOnly, int oneIndexed) const
1350 {
1351 int count = 0;
1352 for(int row=0; row<numRows; row++)
1353 {
1354 if (pointerB != NULL)
1355 pointerB[row] = count + oneIndexed;
1356
1357 int rowLength = GetRowLength(row);
1358 for(int j=0; j< rowLength; j++)
1359 {
1360 if ((!upperTriangleOnly) || (columnIndices[row][j] >= row))
1361 {
1362 if (values != NULL)
1363 values[count] = columnEntries[row][j];
1364 if (columns != NULL)
1365 columns[count] = columnIndices[row][j] + oneIndexed;
1366 count++;
1367 }
1368 }
1369
1370 if (pointerE != NULL)
1371 pointerE[row] = count + oneIndexed;
1372 }
1373 }
1374
Save(const char * filename,int oneIndexed) const1375 int SparseMatrix::Save(const char * filename, int oneIndexed) const
1376 {
1377 FILE * fout = fopen(filename,"w");
1378 if (!fout)
1379 return 1;
1380
1381 fprintf(fout,"%d\n%d\n", numRows, GetNumColumns());
1382 for(int i=0; i<numRows; i++)
1383 {
1384 for(int j=0; j < rowLength[i]; j++)
1385 {
1386 int index = columnIndices[i][j];
1387 double entry = columnEntries[i][j];
1388 fprintf(fout,"%d %d %.15G\n",i + oneIndexed, index + oneIndexed, entry);
1389 }
1390 }
1391 fclose(fout);
1392
1393 return 0;
1394 }
1395
SaveToMatlabFormat(const char * filename) const1396 int SparseMatrix::SaveToMatlabFormat(const char * filename) const
1397 {
1398 FILE * fout = fopen(filename,"w");
1399 if (!fout)
1400 return 1;
1401
1402 for(int i=0; i<numRows; i++)
1403 {
1404 for(int j=0; j < rowLength[i]; j++)
1405 {
1406 int index = columnIndices[i][j];
1407 double entry = columnEntries[i][j];
1408 fprintf(fout,"%d %d %.15G\n",i + 1, index + 1, entry);
1409 }
1410 }
1411 fclose(fout);
1412
1413 return 0;
1414 }
1415
RemoveRowColumn(int index)1416 void SparseMatrix::RemoveRowColumn(int index)
1417 {
1418 // remove row 'index'
1419 free(columnEntries[index]);
1420 free(columnIndices[index]);
1421
1422 for(int i=index; i<numRows-1; i++)
1423 {
1424 columnEntries[i] = columnEntries[i+1];
1425 columnIndices[i] = columnIndices[i+1];
1426 rowLength[i] = rowLength[i+1];
1427 }
1428
1429 // remove column 'index'
1430 for(int i=0; i<numRows-1; i++)
1431 {
1432 // traverse all rows
1433 for(int j=0; j<rowLength[i]; j++)
1434 {
1435 // seek for entry 'index'
1436 if (columnIndices[i][j] == index) // found
1437 {
1438 // shift all elements ahead one step back
1439 for(int k=j; k<rowLength[i]-1; k++)
1440 {
1441 columnIndices[i][k] = columnIndices[i][k+1];
1442 columnEntries[i][k] = columnEntries[i][k+1];
1443 }
1444 rowLength[i]--;
1445 }
1446 }
1447
1448 // decrease indices for DOFs above index
1449 for(int j=0; j<rowLength[i]; j++)
1450 {
1451 if(columnIndices[i][j] > index)
1452 {
1453 // decrease index
1454 columnIndices[i][j]--;
1455 }
1456 }
1457 }
1458
1459 numRows--;
1460 }
1461
RemoveRowsColumnsSlow(int numRemovedRowsColumns,const int * removedRowsColumns,int oneIndexed)1462 void SparseMatrix::RemoveRowsColumnsSlow(int numRemovedRowsColumns, const int * removedRowsColumns, int oneIndexed)
1463 {
1464 for(int i=0; i<numRemovedRowsColumns; i++)
1465 RemoveRowColumn(removedRowsColumns[i]-i-oneIndexed);
1466 }
1467
RemoveRowsColumns(int numRemovedRowsColumns,const int * removedRowsColumns,int oneIndexed)1468 void SparseMatrix::RemoveRowsColumns(int numRemovedRowsColumns, const int * removedRowsColumns, int oneIndexed)
1469 {
1470 // the removed dofs must be pre-sorted
1471 // build a map from old dofs to new ones
1472 vector<int> oldToNew(numRows);
1473 int dof = 0;
1474 int dofCount = 0;
1475 for(int i=0; i<numRemovedRowsColumns; i++)
1476 {
1477 while (dof < removedRowsColumns[i] - oneIndexed)
1478 {
1479 oldToNew[dof] = dofCount;
1480 dofCount++;
1481 dof++;
1482 }
1483 oldToNew[dof] = -1;
1484 dof++;
1485 }
1486 while (dof < numRows)
1487 {
1488 oldToNew[dof] = dofCount;
1489 dofCount++;
1490 dof++;
1491 }
1492
1493 // now, traverse all rows and renumber the entries
1494 int targetRow = 0;
1495 for(int sourceRow = 0; sourceRow < numRows; sourceRow++)
1496 {
1497 if (oldToNew[sourceRow] == -1)
1498 {
1499 free(columnIndices[sourceRow]);
1500 free(columnEntries[sourceRow]);
1501 continue;
1502 }
1503
1504 int targetIndex = 0;
1505 for(int sourceIndex=0; sourceIndex<rowLength[sourceRow]; sourceIndex++)
1506 {
1507 int oldIndex = columnIndices[sourceRow][sourceIndex];
1508 int newIndex = oldToNew[oldIndex];
1509 if (newIndex == -1)
1510 continue;
1511 columnIndices[sourceRow][targetIndex] = newIndex;
1512 columnEntries[sourceRow][targetIndex] = columnEntries[sourceRow][sourceIndex];
1513 targetIndex++;
1514 }
1515
1516 columnIndices[sourceRow] = (int*) realloc(columnIndices[sourceRow], sizeof(int) * targetIndex);
1517 columnEntries[sourceRow] = (double*) realloc(columnEntries[sourceRow], sizeof(double) * targetIndex);
1518
1519 columnIndices[targetRow] = columnIndices[sourceRow];
1520 columnEntries[targetRow] = columnEntries[sourceRow];
1521 rowLength[targetRow] = targetIndex;
1522 targetRow++;
1523 }
1524
1525 numRows -= numRemovedRowsColumns;
1526 columnEntries = (double**) realloc(columnEntries, sizeof(double*) * numRows);
1527 columnIndices = (int**) realloc(columnIndices, sizeof(int*) * numRows);
1528 rowLength = (int*) realloc(rowLength, sizeof(int) * numRows);
1529 }
1530
RemoveColumn(int index)1531 void SparseMatrix::RemoveColumn(int index)
1532 {
1533 // remove column 'index'
1534 for(int i=0; i<numRows; i++)
1535 {
1536 // traverse all rows
1537 for(int j=0; j<rowLength[i]; j++)
1538 {
1539 // seek for entry 'index'
1540 if (columnIndices[i][j] == index) // found
1541 {
1542 // shift all elements ahead one step back
1543 for(int k=j; k<rowLength[i]-1; k++)
1544 {
1545 columnIndices[i][k] = columnIndices[i][k+1];
1546 columnEntries[i][k] = columnEntries[i][k+1];
1547 }
1548 rowLength[i]--;
1549 }
1550 }
1551
1552 // decrease indices for DOFs above index
1553 for(int j=0; j<rowLength[i]; j++)
1554 {
1555 if(columnIndices[i][j] > index)
1556 {
1557 // decrease index
1558 columnIndices[i][j]--;
1559 }
1560 }
1561 }
1562 }
1563
RemoveColumns(int numRemovedColumns,const int * removedColumns,int oneIndexed)1564 void SparseMatrix::RemoveColumns(int numRemovedColumns, const int * removedColumns, int oneIndexed)
1565 {
1566 // the removed dofs must be pre-sorted
1567 // build a map from old dofs to new ones
1568 int numColumns = GetNumColumns();
1569
1570 // must increase numColumns to accommodate matrices with zero columns on the right
1571 for(int i=0; i<numRemovedColumns; i++)
1572 {
1573 int removedColumn0Indexed = removedColumns[i] - oneIndexed;
1574 int neededNumColumns = removedColumn0Indexed + 1;
1575 if (neededNumColumns > numColumns)
1576 numColumns = neededNumColumns;
1577 }
1578
1579 vector<int> oldToNew(numColumns);
1580 int dof = 0;
1581 int dofCount = 0;
1582 for(int i=0; i<numRemovedColumns; i++)
1583 {
1584 while (dof < removedColumns[i] - oneIndexed)
1585 {
1586 oldToNew[dof] = dofCount;
1587 dofCount++;
1588 dof++;
1589 }
1590 oldToNew[dof] = -1;
1591 dof++;
1592 }
1593 while (dof < numColumns)
1594 {
1595 oldToNew[dof] = dofCount;
1596 dofCount++;
1597 dof++;
1598 }
1599
1600 // now, traverse all rows and renumber the entries
1601 for(int row = 0; row < numRows; row++)
1602 {
1603 int targetIndex = 0;
1604 for(int sourceIndex=0; sourceIndex<rowLength[row]; sourceIndex++)
1605 {
1606 int oldIndex = columnIndices[row][sourceIndex];
1607 int newIndex = oldToNew[oldIndex];
1608 if (newIndex == -1)
1609 continue;
1610 columnIndices[row][targetIndex] = newIndex;
1611 columnEntries[row][targetIndex] = columnEntries[row][sourceIndex];
1612 targetIndex++;
1613 }
1614
1615 columnIndices[row] = (int*) realloc(columnIndices[row], sizeof(int) * targetIndex);
1616 columnEntries[row] = (double*) realloc(columnEntries[row], sizeof(double) * targetIndex);
1617 rowLength[row] = targetIndex;
1618 }
1619 }
1620
RemoveColumnsSlow(int numColumns,const int * columns,int oneIndexed)1621 void SparseMatrix::RemoveColumnsSlow(int numColumns, const int * columns, int oneIndexed)
1622 {
1623 for(int i=0; i<numColumns; i++)
1624 RemoveColumn(columns[i]-i-oneIndexed);
1625 }
1626
RemoveRow(int index)1627 void SparseMatrix::RemoveRow(int index)
1628 {
1629 // remove row 'index'
1630 free(columnEntries[index]);
1631 free(columnIndices[index]);
1632
1633 for(int i=index; i<numRows-1; i++)
1634 {
1635 columnEntries[i] = columnEntries[i+1];
1636 columnIndices[i] = columnIndices[i+1];
1637 rowLength[i] = rowLength[i+1];
1638 }
1639
1640 numRows--;
1641 }
1642
RemoveRowsSlow(int numRows,const int * rows,int oneIndexed)1643 void SparseMatrix::RemoveRowsSlow(int numRows, const int * rows, int oneIndexed)
1644 {
1645 for(int i=0; i<numRows; i++)
1646 RemoveRow(rows[i]-i-oneIndexed);
1647 }
1648
RemoveRows(int numRemovedRows,const int * removedRows,int oneIndexed)1649 void SparseMatrix::RemoveRows(int numRemovedRows, const int * removedRows, int oneIndexed)
1650 {
1651 // the removed dofs must be pre-sorted
1652 // build a map from old dofs to new ones
1653 vector<int> oldToNew(numRows);
1654 int dof = 0;
1655 int dofCount = 0;
1656 for(int i=0; i<numRemovedRows; i++)
1657 {
1658 while (dof < removedRows[i] - oneIndexed)
1659 {
1660 oldToNew[dof] = dofCount;
1661 dofCount++;
1662 dof++;
1663 }
1664 oldToNew[dof] = -1;
1665 dof++;
1666 }
1667 while (dof < numRows)
1668 {
1669 oldToNew[dof] = dofCount;
1670 dofCount++;
1671 dof++;
1672 }
1673
1674 // now, traverse all rows and renumber the entries
1675 int targetRow = 0;
1676 for(int sourceRow = 0; sourceRow < numRows; sourceRow++)
1677 {
1678 if (oldToNew[sourceRow] == -1)
1679 {
1680 free(columnIndices[sourceRow]);
1681 free(columnEntries[sourceRow]);
1682 continue;
1683 }
1684
1685 columnIndices[targetRow] = columnIndices[sourceRow];
1686 columnEntries[targetRow] = columnEntries[sourceRow];
1687 rowLength[targetRow] = rowLength[sourceRow];
1688 targetRow++;
1689 }
1690
1691 numRows -= numRemovedRows;
1692 columnEntries = (double**) realloc(columnEntries, sizeof(double*) * numRows);
1693 columnIndices = (int**) realloc(columnIndices, sizeof(double*) * numRows);
1694 rowLength = (int*) realloc(rowLength, sizeof(int) * numRows);
1695 }
1696
GetInfinityNorm() const1697 double SparseMatrix::GetInfinityNorm() const
1698 {
1699 double norm = 0.0;
1700 for(int i=0; i<numRows; i++)
1701 {
1702 double absRowSum = 0;
1703 for(int j=0; j<rowLength[i]; j++)
1704 absRowSum += fabs(columnEntries[i][j]);
1705 if (absRowSum > norm)
1706 norm = absRowSum;
1707 }
1708 return norm;
1709 }
1710
GetFrobeniusNorm() const1711 double SparseMatrix::GetFrobeniusNorm() const
1712 {
1713 double sum = 0.0;
1714 for(int i=0; i<numRows; i++)
1715 {
1716 double rowSum = 0.;
1717 for(int j=0; j<rowLength[i]; j++)
1718 rowSum += columnEntries[i][j] * columnEntries[i][j];
1719 sum += rowSum;
1720 }
1721 return sqrt(sum);
1722 }
1723
DoOneGaussSeidelIteration(double * x,const double * b) const1724 void SparseMatrix::DoOneGaussSeidelIteration(double * x, const double * b) const
1725 {
1726 for(int i=0; i<numRows; i++)
1727 {
1728 double buffer = b[i];
1729 int diagIndex = -1;
1730 for(int j=0; j<rowLength[i]; j++)
1731 {
1732 int column = columnIndices[i][j];
1733 if (column != i)
1734 buffer -= columnEntries[i][j] * x[column];
1735 else
1736 diagIndex = j;
1737 }
1738 x[i] = buffer / columnEntries[i][diagIndex];
1739 }
1740 }
1741
ComputeResidual(const double * x,const double * b,double * residual) const1742 void SparseMatrix::ComputeResidual(const double * x, const double * b, double * residual) const
1743 {
1744 MultiplyVector(x,residual);
1745 for(int i=0; i<numRows; i++)
1746 residual[i] -= b[i];
1747 }
1748
CheckLinearSystemSolution(const double * x,const double * b,int verbose,double * buffer) const1749 double SparseMatrix::CheckLinearSystemSolution(const double * x, const double * b, int verbose, double * buffer) const
1750 {
1751 double * bufferv = NULL;
1752
1753 if (buffer == NULL)
1754 {
1755 bufferv = (double*) malloc (sizeof(double) * numRows);
1756 buffer = bufferv;
1757 }
1758
1759 MultiplyVector(x,buffer);
1760
1761 double inftyNorm = 0;
1762 double inftyNorm_b = 0;
1763 for(int i=0; i<numRows; i++)
1764 {
1765 if (fabs(buffer[i] - b[i]) > inftyNorm)
1766 inftyNorm = fabs(buffer[i] - b[i]);
1767
1768 if (fabs(b[i]) > inftyNorm_b)
1769 inftyNorm_b = fabs(b[i]);
1770 }
1771
1772 if (verbose)
1773 {
1774 printf("Infinity residual norm ||Ax-b|| is %G. ||b|| is %G.\n", inftyNorm, inftyNorm_b);
1775 printf("Relative infinity residual norm ||Ax-b||/||b|| is %G.\n", inftyNorm / inftyNorm_b);
1776 }
1777
1778 free(bufferv);
1779
1780 return inftyNorm / inftyNorm_b;
1781 }
1782
MakeDenseMatrix(double * denseMatrix,int numColumns) const1783 void SparseMatrix::MakeDenseMatrix(double * denseMatrix, int numColumns) const
1784 {
1785 if (numColumns == -1)
1786 numColumns = GetNumColumns();
1787 memset(denseMatrix, 0, sizeof(double) * (numRows * numColumns));
1788 for(int i=0; i< numRows; i++)
1789 for(int j=0; j<rowLength[i]; j++)
1790 {
1791 if (columnIndices[i][j] < numColumns)
1792 denseMatrix[numRows * columnIndices[i][j] + i] = columnEntries[i][j];
1793 }
1794 }
1795
MakeDenseMatrixTranspose(int numColumns,double * denseMatrix) const1796 void SparseMatrix::MakeDenseMatrixTranspose(int numColumns, double * denseMatrix) const
1797 {
1798 // note: we cannot use GetNumColumns() here because the rightmost columns of the sparse matrix can be zero and the GetNumColumns() will not be accurate
1799 memset(denseMatrix, 0, sizeof(double) * (numRows * numColumns));
1800 for(int i=0; i<numRows; i++)
1801 {
1802 int offset = i * numColumns;
1803 for(int j=0; j<rowLength[i]; j++)
1804 denseMatrix[offset + columnIndices[i][j]] = columnEntries[i][j];
1805 }
1806 }
1807
MultiplyRow(int row,double scalar)1808 void SparseMatrix::MultiplyRow(int row, double scalar) // multiplies all elements in row 'row' with scalar 'scalar'
1809 {
1810 for(int j=0; j<rowLength[row]; j++)
1811 columnEntries[row][j] *= scalar;
1812 }
1813
GetNumColumns() const1814 int SparseMatrix::GetNumColumns() const
1815 {
1816 int numColumns = -1;
1817 for(int i=0; i<numRows; i++)
1818 {
1819 for(int j=0; j<rowLength[i]; j++)
1820 {
1821 if (columnIndices[i][j] > numColumns)
1822 numColumns = columnIndices[i][j];
1823 }
1824 }
1825 return numColumns + 1;
1826 }
1827
IncreaseNumRows(int numAddedRows)1828 void SparseMatrix::IncreaseNumRows(int numAddedRows)
1829 {
1830 int newn = numRows + numAddedRows;
1831
1832 rowLength = (int*) realloc (rowLength, sizeof(int) * newn);
1833 for(int i=0; i<numAddedRows; i++)
1834 rowLength[numRows + i] = 0;
1835
1836 columnIndices = (int**) realloc (columnIndices, sizeof(int*) * newn);
1837 for(int i=0; i<numAddedRows; i++)
1838 columnIndices[numRows + i] = NULL;
1839
1840 columnEntries = (double**) realloc (columnEntries, sizeof(double*) * newn);
1841 for(int i=0; i<numAddedRows; i++)
1842 columnEntries[numRows + i] = NULL;
1843
1844 numRows = newn;
1845 }
1846
ConjugateMatrix(SparseMatrix & U,int verbose,int numColumns)1847 SparseMatrix SparseMatrix::ConjugateMatrix(SparseMatrix & U, int verbose, int numColumns)
1848 {
1849 if (numColumns == -1)
1850 numColumns = U.GetNumColumns();
1851
1852 SparseMatrixOutline outline(numColumns);
1853
1854 for(int i=0; i<numRows; i++)
1855 {
1856 if (verbose)
1857 {
1858 if (i % 100 == 1)
1859 printf("Processing row %d / %d...\n", i, numRows);
1860 }
1861 for(int j=0; j<rowLength[i]; j++)
1862 {
1863 int I = i;
1864 int J = columnIndices[i][j];
1865 double scalar = columnEntries[i][j];
1866
1867 // compute tensor product of rows I and J of U
1868 for(int k=0; k<U.rowLength[I]; k++)
1869 for(int l=0; l<U.rowLength[J]; l++)
1870 {
1871 int K = U.columnIndices[I][k];
1872 int L = U.columnIndices[J][l];
1873
1874 // there is an entry at (I, K), and another entry at (J, L); compute their contribution to tensor product:
1875 outline.AddEntry(K, L, scalar * U.columnEntries[I][k] * U.columnEntries[J][l]);
1876 }
1877 }
1878 }
1879
1880 if (verbose)
1881 printf("Creating sparse matrix from outline...\n");
1882
1883 return SparseMatrix(&outline);
1884 }
1885
BuildConjugationIndices(SparseMatrix & U,SparseMatrix & MTilde,precomputedIndicesType * precomputedIndices)1886 void SparseMatrix::BuildConjugationIndices(SparseMatrix & U, SparseMatrix & MTilde, precomputedIndicesType * precomputedIndices)
1887 {
1888 typedef pair< pair<int,int>, pair<int, int> > fourTuple;
1889 typedef vector<fourTuple> listOfFourTuples;
1890 typedef map<int, listOfFourTuples> rowMap;
1891 vector<rowMap> rowMaps;
1892 for(int i=0; i<MTilde.numRows; i++)
1893 rowMaps.push_back(rowMap());
1894
1895 for(int i=0; i<numRows; i++)
1896 {
1897 for(int j=0; j<rowLength[i]; j++)
1898 {
1899 int I = i;
1900 int J = columnIndices[i][j];
1901
1902 // compute tensor product of rows I and J of U
1903 for(int k=0; k<U.rowLength[I]; k++)
1904 for(int l=0; l<U.rowLength[J]; l++)
1905 {
1906 int K = U.columnIndices[I][k];
1907 int L = U.columnIndices[J][l];
1908 //outline.AddEntry(K, L, scalar * U.columnEntries[I][k] * U.columnEntries[J][l]);
1909
1910 fourTuple tuple(make_pair(make_pair(i,j), make_pair(k,l)));
1911
1912 rowMap::iterator iter = rowMaps[K].find(L);
1913 if (iter == rowMaps[K].end())
1914 {
1915 listOfFourTuples singletonList;
1916 singletonList.push_back(tuple);
1917 rowMaps[K].insert(make_pair(L, singletonList));
1918 }
1919 else
1920 {
1921 (iter->second).push_back(tuple);
1922 }
1923 }
1924 }
1925 }
1926
1927 // copy map to precomputedIndices
1928 (*precomputedIndices) = (int***) malloc (sizeof(int**) * MTilde.numRows);
1929 for(int i=0; i<MTilde.numRows; i++)
1930 {
1931 (*precomputedIndices)[i] = (int**) malloc (sizeof(int*) * rowMaps[i].size());
1932 int j = 0;
1933 for(rowMap::iterator iter = rowMaps[i].begin(); iter != rowMaps[i].end(); ++iter)
1934 {
1935 (*precomputedIndices)[i][j] = (int*) malloc (sizeof(int) * (4 * ((iter->second).size()) + 1));
1936 ((*precomputedIndices)[i][j])[0] = (int)((iter->second).size());
1937 for(int k=0; k<((*precomputedIndices)[i][j])[0]; k++)
1938 {
1939 ((*precomputedIndices)[i][j])[1+4*k+0] = ((iter->second)[k]).first.first;
1940 ((*precomputedIndices)[i][j])[1+4*k+1] = ((iter->second)[k]).first.second;
1941 ((*precomputedIndices)[i][j])[1+4*k+2] = ((iter->second)[k]).second.first;
1942 ((*precomputedIndices)[i][j])[1+4*k+3] = ((iter->second)[k]).second.second;
1943 }
1944 j++;
1945 }
1946 }
1947 }
1948
ConjugateMatrix(precomputedIndicesType precomputedIndices,SparseMatrix & U,SparseMatrix & MTilde)1949 void SparseMatrix::ConjugateMatrix(precomputedIndicesType precomputedIndices, SparseMatrix & U, SparseMatrix & MTilde)
1950 {
1951 MTilde.ResetToZero();
1952 for(int row=0; row<MTilde.numRows; row++)
1953 {
1954 int ** rowIndices = precomputedIndices[row];
1955 for(int j=0; j<MTilde.rowLength[row]; j++)
1956 {
1957 int * entryIndices = rowIndices[j];
1958 int numSummationTerms = entryIndices[0];
1959 for(int k=0; k<numSummationTerms; k++)
1960 {
1961 int * entryIndex = &entryIndices[4*k+1];
1962 int rowOfM = entryIndex[0];
1963 int columnIndexOfM = entryIndex[1];
1964 int columnOfM = columnIndices[rowOfM][columnIndexOfM];
1965 int columnIndexofU_for_MTilde_row = entryIndex[2];
1966 int columnIndexofU_for_MTilde_column = entryIndex[3];
1967 (MTilde.columnEntries)[row][j] += columnEntries[rowOfM][columnIndexOfM] * U.columnEntries[row][columnIndexofU_for_MTilde_row] * U.columnEntries[columnOfM][columnIndexofU_for_MTilde_column];
1968 }
1969 }
1970 }
1971 }
1972
ConjugateMatrix(const double * U,int r,double * UTilde) const1973 void SparseMatrix::ConjugateMatrix(const double * U, int r, double * UTilde) const
1974 {
1975 double * MU = (double*) malloc (sizeof(double) * numRows * r);
1976 MultiplyMatrix(numRows, r, U, MU);
1977
1978 // compute U^T * MU
1979 for(int i=0; i<r; i++)
1980 for(int j=0; j<r; j++)
1981 {
1982 double entry = 0.0;
1983 for(int k=0; k<numRows; k++)
1984 entry += U[i * numRows + k] * MU[j * numRows + k];
1985 UTilde[j * r + i] = entry;
1986 }
1987
1988 free(MU);
1989 }
1990
Transpose(int numColumns)1991 SparseMatrix * SparseMatrix::Transpose(int numColumns)
1992 {
1993 if (numColumns < 0)
1994 numColumns = GetNumColumns();
1995
1996 SparseMatrixOutline outline(numColumns);
1997
1998 for(int i=0; i<numRows; i++)
1999 for(int j=0; j<rowLength[i]; j++)
2000 outline.AddEntry(columnIndices[i][j], i, columnEntries[i][j]);
2001
2002 return new SparseMatrix(&outline);
2003 }
2004
AssignTransposedMatrix(SparseMatrix & AT)2005 void SparseMatrix::AssignTransposedMatrix(SparseMatrix & AT)
2006 {
2007 AT.BuildTranspositionIndices();
2008 for(int i = 0; i < AT.numRows; i++)
2009 for(int j = 0; j < AT.rowLength[i]; j++)
2010 {
2011 int index = AT.TransposedIndex(i,j);
2012 int row = AT.columnIndices[i][j];
2013 columnEntries[row][index] = AT.columnEntries[i][j];
2014 }
2015 }
2016
SetRows(const SparseMatrix * source,int startRow,int startColumn)2017 void SparseMatrix::SetRows(const SparseMatrix * source, int startRow, int startColumn)
2018 {
2019 for(int i=0; i<source->GetNumRows(); i++)
2020 {
2021 int row = startRow + i;
2022 if (row >= numRows)
2023 return;
2024
2025 rowLength[row] = source->GetRowLength(i);
2026 columnIndices[row] = (int*) realloc (columnIndices[row], sizeof(int) * rowLength[row]);
2027 columnEntries[row] = (double*) realloc (columnEntries[row], sizeof(double) * rowLength[row]);
2028 for(int j=0; j<rowLength[row]; j++)
2029 {
2030 columnIndices[row][j] = startColumn + source->columnIndices[i][j];
2031 columnEntries[row][j] = source->columnEntries[i][j];
2032 }
2033 }
2034 }
2035
AppendRowsColumns(const SparseMatrix * source)2036 void SparseMatrix::AppendRowsColumns(const SparseMatrix * source)
2037 {
2038 int * oldRowLengths = (int*) malloc (sizeof(int) * numRows);
2039 for(int i=0; i<numRows; i++)
2040 oldRowLengths[i] = rowLength[i];
2041
2042 int oldNumRows = numRows;
2043 IncreaseNumRows(source->GetNumRows());
2044 SetRows(source, oldNumRows);
2045
2046 // add transpose of rows:
2047
2048 // first, establish new column lengths
2049 for(int row=0; row<source->GetNumRows(); row++)
2050 {
2051 for(int j=0; j<source->GetRowLength(row); j++)
2052 {
2053 int column = source->GetColumnIndex(row, j);
2054 rowLength[column]++;
2055 }
2056 }
2057
2058 // extend size
2059 for(int row=0; row<oldNumRows; row++)
2060 {
2061 columnIndices[row] = (int*) realloc (columnIndices[row], sizeof(int) * rowLength[row]);
2062 columnEntries[row] = (double*) realloc (columnEntries[row], sizeof(double) * rowLength[row]);
2063 }
2064
2065 // write entries into their place
2066 for(int i=0; i<oldNumRows; i++)
2067 rowLength[i] = oldRowLengths[i];
2068
2069 for(int row=0; row<source->GetNumRows(); row++)
2070 {
2071 for(int j=0; j<source->GetRowLength(row); j++)
2072 {
2073 int column = source->GetColumnIndex(row, j);
2074 columnIndices[column][rowLength[column]] = oldNumRows + row;
2075 columnEntries[column][rowLength[column]] = source->GetEntry(row, j);
2076 rowLength[column]++;
2077 }
2078 }
2079
2080 free(oldRowLengths);
2081
2082 // append zero diagonal in lower-right block (helps with some solvers, such as PARDISO)
2083 for(int row=0; row<source->GetNumRows(); row++)
2084 {
2085 rowLength[oldNumRows + row]++;
2086 columnIndices[oldNumRows + row] = (int*) realloc (columnIndices[oldNumRows + row], sizeof(int) * rowLength[oldNumRows + row]);
2087 columnEntries[oldNumRows + row] = (double*) realloc (columnEntries[oldNumRows + row], sizeof(double) * rowLength[oldNumRows + row]);
2088 columnIndices[oldNumRows + row][rowLength[oldNumRows + row] - 1] = oldNumRows + row;
2089 columnEntries[oldNumRows + row][rowLength[oldNumRows + row] - 1] = 0.0;
2090 }
2091 }
2092
CreateIdentityMatrix(int numRows)2093 SparseMatrix * SparseMatrix::CreateIdentityMatrix(int numRows)
2094 {
2095 SparseMatrixOutline * outline = new SparseMatrixOutline(numRows);
2096 for (int row=0; row<numRows; row++)
2097 outline->AddEntry(row, row, 1.0);
2098 SparseMatrix * mat = new SparseMatrix(outline);
2099 delete outline;
2100 return mat;
2101 }
2102