C++ Help Required

[sibarrick]

Hi,

I'm trying to use UT_SparseMatrix http://odforce.net/hdk/tree/H8.0/toolkit/html/ but I can't see how you access an individual cell value once you have set up the matrix. I'm thinking maybe it has something to do with the ut_MatrixCell class that is defined within the UT_SparseMatrix class, but this is beyond my C++ knowledge at present. Can a C++ guru explain how you do this?

Also rather worryingly the header for H9 shows this class has totally changed..... so even if I get this to work in H8 it looks like it will then fail in H9

[edward]

The quick answer is you don't. It's all geared towards matrix solving right now. If you must, you can use extractSubMatrix() and then multiply the result by a UT_Vector of 1. However, it's likely gonna be slow.

[Mario Marengo]

Yeah, it looks like it was built for quickly applying operators, not for randomly accessing/editing elements. On the other hand, if you really want to live dangerously, break all the rules, and guarantee yourself a choice location in "C++ Hell", you could always hack your way in like this:

#include <iostream>
#include <UT/UT_SparseMatrix.h>

struct hackdata
{
   int	  myRow;
   int	  myCol;
   fpreal64 myValue;
};

struct hack 
{
   int			   myNumRows;
   int			   myNumCols;
   mutable hackdata *myCells;
   mutable int	   myCount;
   int			   myMaxSize;
   mutable bool	  myCompiledFlag;
};


// Using the definition in the H9 headers...
int main() {
   UT_SparseMatrix M(2,2);
   M.addToElement(0,0,1); // Create a non-zero entry
   if(M.getNumRows()) {
	  hackdata *cells = reinterpret_cast<hack *>(&M)->myCells;
	  std::cout<<"Data at row "<<cells[0].myRow
			   <<", column "<<cells[0].myCol
			   <<", is "<<cells[0].myValue<<std::endl;
   }
   return 0;
}

The above just clones the static data of the class and plays with pointers to gain access... and it's all kinds of bad for your health. But fun to play with (note it uses the class definitions in the H9 headers).

Be sure to send us a postcard from the Dark Side!

Cheers.

[sibarrick]

Hmmmm.... not wishing to cross to the dark side I think I'll have to take a different path.

I think the best thing for me to do is store the values I need that I'm currently poking into the matrix in a seperate vector at the same time. I only need to access them once. The way the code is arranged it would be easier to access the values I need afterwards, but I'm sure I can find a way to stuff the ones I need into a vector as I go rather than get too hacky. Even I have my limits

[peship]

This one works really well: http://math.nist.gov/sparselib++/

[sibarrick]

cheers. I'm going to try and do it with the built in hdk solvers first as I know a bit more what I'm doing with them.

A related question though. If you are building a diagonal (jacobi??) matrix pre-conditioner from a non-symmetric matrix, what constitutes the diagonal??

| 1 0 0 0 0 0 |

| 0 1 0 0 0 0 |

| 0 0 1 0 0 0 |

| 0 0 0 1 0 0 |

| 0 0 0 0 1 0 |

ie the i, i diagonal or

something complicated that somehow takes into account the non-symmetry so the diagonal reaches the top-left to bottom-right positions?

| 1 0 0 0 0 0 |

| 0 1 0 0 0 0 |

| 0 0 1 0 0 0 |

| 0 0 0 0 1 0 |

| 0 0 0 0 0 1 |

or does it actually not really matter since its just a starting point for the solver anyhow??

[Mario Marengo]

A related question though. If you are building a diagonal (jacobi??) matrix pre-conditioner from a non-symmetric matrix, what constitutes the diagonal??

Hey Simon,

I don't know the answer to your question(s), but I was under the impression that the property of symmetry (or non-symmetry), as well as diagonalization was only applicable to square matrices. i.e: a matrix can be non-symmetric in the numerical sense and still be square. For example, here is a symmetric matrix (it's identical to its transpose):

|a b|

|b a|

And a non-symmetric one (it's different than its transpose):

|a b|

|c a|

Yet they are both square.

For the rest of it, I don't know what your context is, but as I understand it, the "preconditioner" of a linear system (matrix) is another matrix that is used to modify the original so that, when some iterative method is used to solve the system, it will converge on a solution faster. The "Jacobi Preconditioner" is one method for arriving at this preconditioning matrix, and it boils down to simply keeping the diagonal elements of the original matrix (i.e: setting everything else to 0), and taking its inverse (which in the case of a diagonal matrix is simply the numerical inverse of each one of its diagonal elements).

So for a matrix like this one (which is non-symmetric, and not decomposed into any canonical form):

|a b|

|c d|

The Jacobi Preconitioner matrix would be:

|1/a 0|

|0 1/d|

However, if in doubt, I'm sure you could fall back to the identity matrix as a preconditioner, although this would be the same as trying to solve the original system "without any help", if you see what I mean.

Hopefully someone who really knows this stuff can set us both straight though

Cheers.

[sibarrick]

Sorry for the confusion, i keep getting my matrix terminology mixed up, I'll try and elaborate.

I was talking about a matrix that is non-square and non-symmetric.

Where I put the 1's i probably should have used a,b,c's i didn't mean to imply it was an identy matrix, i was just indicating positions with the 1's.

As far as I understand it (which admittedly so far isn't very far) there are methods for solving non-square/non-symmetric matrices, there is one in the hdk PCGLS which requires a preconditioner to work. Again as I understand it you can supply a Jacobi Preconditioner to this method. So what I should have said is what is the definition of a diagonal in a non-square matrix, not just a non-symmetric matrix so I can set up such a preconditioner.

Always hated matrices (never meet anyone who didn't) but man they are powerful solvers..... just not design for human brains.

[edward]

So why not do what Mario suggests and simply not precondition your matrix? Your callback can just assign b to x. But I think it's the first case except it assumes that x is of size n (not m).

[sibarrick]

Well I'll definetely try that, I was just rather hoping there was a way to do it.

[edward]

Sorry, the AtASolve() is actually solving for A^T.A (A transpose times A). So it should be preconditioning the square matrix, A^T.A which is of size n by n. So eg.

static void
AtASolve(const UT_Vector &b, UT_Vector &x)
{
	for (int i = 0; i <= x.getNH(); i++)
	{
		if (AtADiag(i)) // check for divide by 0
			x(i) = b(i) / AtADiag(i);
		else
			x(i) = b(i);
	}
}

[sibarrick]

I see. I guess that brings me back to the beginning though, once i have A^T.A I then need to extract the diagonal members which as we have established is rather tricky using UT_SparseMatrix. I guess I could do that bit in a non sparse matrix form but that would kind of defeat the object.... maybe I'll start off that way just to see if I can get it working and then worry about efficiency later.

[edward]

Ah, but UT_SparseMatrix has a method to compute this directly ...

AtADiag.init(0, A.getNumCols()-1);
A.allColNorm2(AtADiag);

[sibarrick]

Ah hah! it all drops into place!!!!