[jasper_db]

Hi,

I am a student civil engineering. I am using Xcode to create a little tool. In this tool I need to solve sets of equations and use other matrix operations. I implemented matrix multiplications and a simple Gauss-eliminations with backsubstitution for solving the equations as objective-c code without the help of other frameworks.

Now after a search on the web I found that the accelerate-framework from apple deals with these things. And that it implements the Lapack-routines. After a little search a don't seem to find how I can use this framework. Is there someone who can give me a simple code-example for using this framework to solve a set of equations A*x=B?

greetings
Jasper

[MarkLT1]

jasper-

What type of Civ E work are you doing?

The use of blas/lapack in c/c++/obj c can be a tricky endeavor, as LAPACK/BLAS were originally written for fortran. That being said, Apple has made it somewhat easy by also including a cblas.h and clapack.h header file in the accelerate framework. The following code is a simple code, using dynamically allocated c++ arrays. You most likely would want to use some sort of vector/matrix class (or in c++ use the stl vector class as your storage scheme). A few gotchas-

1) Fortran assumes you are using a major column ordered matrix, and thus blas/lapack assume the same. So your data structure would look like:
{row1_col1, row2_col1, ... , rowM_col1, row1_col2, row2_col2, ... , rowM_col2, ... row1_colN, row2_colN, ... ,rowM_colN}

2) if you get into using sparse storage schemes, Fortran assumes the first index is 1, while C assumes it is zero.. so you need to keep this in mind.

So here is the sample C++ code (hopefully this will look OK). It would be compiled at the command line via:

g++ file.cpp -o file -framework accelerate

[CODEX]
#include <math.h>
#include <iostream>

// include the clapack headers
#include <vecLib/clapack.h>

using namespace std;

int main (int argc, char const *argv[])
{
// Sample problem to solve Ax=b for x. System of equations will be:
// A=[[1. 2. 3. 4.]
// [5. 6. 7. 8.]
// [9. 10. 11. 12.]
// [13. 14. 15. 16.]]
//
// b = [17.
// 18.
// 19.
// 20.]
// Dimensions
int num_rows=4;
int num_cols=4;

// create data array for our matrix. Will be a vector that is num_rows*num_cols elements
// to access A(r,c) use A[r+c*num_rows], in otherwords major column ordering. {r1c1,r2c1,r3c1...}
double * A;
A = new double[num_rows*num_cols];

// Zero out A;
for(int i=0; i<num_rows*num_cols; i++)
A[i]=0.0;

// Fill in A the "normal" way
double val=1.0;
for(int row=0; row<num_rows; row++) {
for(int col=0; col<num_cols; col++) {
A[row+col*num_rows]=val;
val+=1.0;
}
}

// display A:
cout << "A=" << endl;
for(int row=0; row<num_rows; row++) {
for(int col=0; col<num_cols; col++) {
cout << A[row+col*num_rows] << " ";
}
cout << endl;
}
cout << endl;

// Create data array for b. Will be a vector that is num_rows in size.
double *b;
b = new double[num_rows];

// zero out B:
for(int i=0; i<num_rows; i++)
b[i]=0.0;

// fill in B:
for(int i=0; i<num_rows; i++) {
b[i]=val;
val+=1.0;
}

// display b.
cout << "b=" << endl;
for(int i=0; i<num_rows; i++) {
cout << b[i] << endl;
}
cout << endl;

// create x vector with which will be our solution vector
double *x = new double[num_rows];
for(int i=0; i<num_rows; i++) {
x[i]=b[i]; // copy b into x, as the lapack routines do an inplace solve
} // i.e.- the rhs vector input, is overwritten with the solution vector
// so if we want to preserve b, we need to make a copy.

//-------------------------//
// The Solver- using dgesv //
//-------------------------//
__CLPK_integer n, lda, ldb, nrhs, info; // Need some parameters/storage vars of type clapack integer
n=lda=ldb=num_rows; // set n, lda, ldb to num_rows
nrhs=1; // we are solving for 1 rhs vector
__CLPK_integer * ipiv = new __CLPK_integer[num_rows]; // create a pivot vector for the solver
dgesv_(&n, &nrhs, A, &lda, ipiv, x, &ldb, &info); // perform the [d]ouble precision [ge]eneral Matrix [s]ol[v]e (dgesv)

// output x
cout << "x=" << endl;
for(int i=0; i<num_rows; i++) {
cout << x[i] << endl;
}
return 0;
}
[/CODEX]

[jasper_db]

First of all, thank you for helping me out! I am going to anayse the code a little. But it seems to work fine at first sight.

Would I encounter problems because of the c++ if I woud try to integrate this in my current objective-c foundation tool?

I needed this for the following reason. Last semester I had a course about the finite element-method. In the little time I have in between, I am trying to build my own simple (2D, only beam elements) FE-application for the structural static analysis of a building. So this means a system of equations KU=P needs to be solved.

Maybe the use of Lapack for a newbie is quite ambitious... But otherwise I have to implement the algorithms myself.

thx again
Jasper

[MarkLT1]

Originally Posted by jasper_db 

First of all, thank you for helping me out! I am going to anayse the code a little. But it seems to work fine at first sight.

No problem. I hope it helps!

Would I encounter problems because of the c++ if I woud try to integrate this in my current objective-c foundation tool?

I have never programmed in Obj. C myself, but seeing as it is derived from (and is supposed to be very similar to) C, I would guess the code wouldn't change too much. The basic premise is that the Lapack routine I used takes as input some integers that describe the data, and 3 vectors- a vector of doubles that describes the coefficient matrix A (in major column ordering), a vector of doubles that describes the right hand side, b (and will also be replaced with the solution vector), and a vector of integers that is used as a pivot table.

I needed this for the following reason. Last semester I had a course about the finite element-method. In the little time I have in between, I am trying to build my own simple (2D, only beam elements) FE-application for the structural static analysis of a building. So this means a system of equations KU=P needs to be solved.

Very cool. That is the type of work I do all day. For my dissertation work, I was doing non-linear dynamic time-history modeling of structures subjected to earthquakes. I have now moved over to the environmental engineering side, and do a bunch of flow through porous media using FEM and Finite Volume methods.

Maybe the use of Lapack for a newbie is quite ambitious... But otherwise I have to implement the algorithms myself.

Nah.. you should be fine. The lapack routines are just function calls- as long as you know the data that needs to be fed in, and how it is handled, it is cake. Take a look through the clapack.h header file in the Accelerate Framework for more detailed info on the functions and inputs.

[Chuckit]

Objective-C on the Mac is pretty much fully compatible with C++.

[jasper_db]

Thx all for the reactions! Really cool to have some guys who want help me!

As the documentation said, just change the extension of the file from .cpp to .mm, add the framework to my project and I can add objective-c to your code and everything compiles well in Xcode.

Funny you have worked on dynamic loadings on structures. There is a big chance I will spend my thesis next year on that subject! :-)

thx again for the help!
Jasper

[MarkLT1]

Originally Posted by jasper_db 

Funny you have worked on dynamic loadings on structures. There is a big chance I will spend my thesis next year on that subject! :-)

In that case, dont re-invent the wheel. Check out OpenSees (Open System for Earthquake Engineering Simulation - Home Page). It is a command line Dynamic Structural Analysis program, and is AMAZING. Compiles pretty easily on OS X.

[jasper_db]

Just to say I integrated the Lapack-solver successful in my app.

As I am not familiar with C++. The way to deallocate the vector A and b is just to put "delete A; delete b;", not?

[MarkLT1]

jasper-

Congrats on getting your code to work!

In C++, to deallocate an array you use the delete function, but also have to indicate an array with []. As you currently have your code (delete a; delete b) you are not specifying that the memory to be freed is an array, so you will have problems. Just insert [] as in the code below. Also, arrays are simply a pointer, so it is a good idea to set the variable to zero (the null pointer) so that the variable is not pointing to a non-existent spot in memory.

[CODEX]
delete [] myArray; // Free the memory associated with myArray
myArray = 0; // Set myArray to the null pointer. Could use NULL in place of 0. This is technically not necessary, but is good, clean, practice.
[/CODEX]

edit: If you are doing anything more complicated than this simple example, you will probably want to use some sort of wrapper around the C++ arrays, by either writing your own (that does its own cleanup) or by using the STL::vector library.

Lost of good info on the STL vector lib here: http://www.cplusplus.com/reference/stl/vector/

For an STL Vector, when you make your call to lapack, you'd use the myvector.front() to reference the first element of your array instead of the variable name. If this is confusing, and you want to know more, let me know, and I'll post some example code later.

[jasper_db]

I ain't all that serious :-).

The basic FE-app is actually just the engine of the GUI-app I wanted to create. Just wanted to learn to build a Cocoa-app on the mac.
But who knows how this little experiment will evolve later...

But anyway thank you very much for all the help! Although you can learn a lot by yourself, but sometimes you just need a little help from the outside!

Jasper