Jump to content
Sign in to follow this  
sibarrick

C++ Help Required

Recommended Posts

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

Share this post


Link to post
Share on other sites

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. :)

Share this post


Link to post
Share on other sites

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.

Share this post


Link to post
Share on other sites

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

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 :)

Share this post


Link to post
Share on other sites

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??

Share this post


Link to post
Share on other sites
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.

Share this post


Link to post
Share on other sites

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.

Share this post


Link to post
Share on other sites

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).

Share this post


Link to post
Share on other sites

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);
	}
}

Share this post


Link to post
Share on other sites

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.

Share this post


Link to post
Share on other sites

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

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

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now
Sign in to follow this  

×