paso/src/SparseMatrix.c


/*******************************************************
*
* Copyright (c) 2003-2008 by University of Queensland
* Earth Systems Science Computational Center (ESSCC)
* http://www.uq.edu.au/esscc
*
* Primary Business: Queensland, Australia
* Licensed under the Open Software License version 3.0
* http://www.opensource.org/licenses/osl-3.0.php
*
*******************************************************/


/**************************************************************/

/* Paso: SparseMatrix */

/**************************************************************/

/* Author: gross@access.edu.au */

/**************************************************************/

#include "Paso.h"
#include "SparseMatrix.h"
#include "TRILINOS.h"
#include "mmio.h"

/**************************************************************/
static void swap( index_t*, index_t*, double*, int, int );
static void q_sort( index_t*, index_t*, double*, int, int );
/*static void print_entries( index_t*, index_t*, double* );*/

static int M, N, nz;


/* debug: print the entries */
/*
void print_entries( index_t *r, index_t *c, double *v )
{
    int i;

    for( i=0; i<nz; i++ )
    {
        printf( "(%ld, %ld) == %e\n", (long)r[i], (long)c[i], v[i] );
    }
}
*/

/* swap function */
void swap( index_t *r, index_t *c, double *v, int left, int right )
{
    double v_temp;
    index_t temp;

    temp = r[left];
    r[left] = r[right];
    r[right] = temp;

    temp = c[left];
    c[left] = c[right];
    c[right] = temp;

    v_temp = v[left];
    v[left] = v[right];
    v[right] = v_temp;
}

void q_sort( index_t *row, index_t *col, double *val, int begin, int end )
{
    int l, r;
    index_t pivot, lval;

    if( end > begin )
    {
        pivot = N * row[begin] + col[begin];
        l = begin + 1;
        r = end;

        while( l < r )
        {
            lval = N * row[l] + col[l];
            if( lval < pivot )
                l++;
            else
            {
                r--;
                swap( row, col, val, l, r );
            }
        }
        l--;
        swap( row, col, val, begin, l );
        q_sort( row, col, val, begin, l );
        q_sort( row, col, val, r, end );
    }
}


/* allocates a SparseMatrix of type type using the given matrix pattern 
   if type is UNKOWN CSR is used.
   if CSC or CSC_BLK1 is used pattern has to give the CSC pattern.
   if CSR or CSR_BLK1 is used pattern has to give the CSR pattern.
   Values are initialized by zero.  */
Paso_SparseMatrix* Paso_SparseMatrix_alloc(Paso_SparseMatrixType type,Paso_Pattern *pattern, int row_block_size, int col_block_size) {

  double time0;
  Paso_SparseMatrix*out=NULL;
  Paso_SparseMatrixType pattern_format_out;

  Paso_resetError();
  time0=Paso_timer();
  out=MEMALLOC(1,Paso_SparseMatrix);
  if (! Paso_checkPtr(out)) {  
     out->pattern=NULL;  
     out->val=NULL;  
     out->reference_counter=1;
     out->type=type;

     pattern_format_out= (type & MATRIX_FORMAT_OFFSET1)? PATTERN_FORMAT_OFFSET1:  PATTERN_FORMAT_DEFAULT;
     /* ====== compressed sparse columns === */
     if (type & MATRIX_FORMAT_CSC) {
        if (type & MATRIX_FORMAT_SYM) {
           Paso_setError(TYPE_ERROR,"Generation of matrix pattern for symmetric CSC is not implemented yet.");
           return NULL;
        } else {
           if ((type & MATRIX_FORMAT_BLK1) || row_block_size!=col_block_size || col_block_size>3) {
              out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,col_block_size,row_block_size);
              out->row_block_size=1;
              out->col_block_size=1;
           } else {
             out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,1,1);
              out->row_block_size=row_block_size;
              out->col_block_size=col_block_size;
           }
           out->numRows = out->pattern->numInput;
           out->numCols = out->pattern->numOutput;
        }
     } else {
     /* ====== compressed sparse row === */
        if (type & MATRIX_FORMAT_SYM) {
           Paso_setError(TYPE_ERROR,"Generation of matrix pattern for symmetric CSR is not implemented yet.");
           return NULL;
        } else {
           if ((type & MATRIX_FORMAT_BLK1) || row_block_size!=col_block_size || col_block_size>3)  {
              out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,row_block_size,col_block_size);
              out->row_block_size=1;
              out->col_block_size=1;
           } else {
              out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,1,1);
              out->row_block_size=row_block_size;
              out->col_block_size=col_block_size;
           }
           out->numRows = out->pattern->numOutput;
           out->numCols = out->pattern->numInput;
        }
     }
     out->block_size=out->row_block_size*out->col_block_size;
     out->len=(size_t)(out->pattern->len)*(size_t)(out->block_size);

     out->val=MEMALLOC(out->len,double);
     if (! Paso_checkPtr(out->val)) Paso_SparseMatrix_setValues(out,DBLE(0));
  }
  /* all done: */
  if (! Paso_noError()) {
    Paso_SparseMatrix_free(out);
    return NULL;
  } else {
    #ifdef Paso_TRACE
    printf("timing: system matrix %.4e sec\n",Paso_timer()-time0);
    printf("Paso_SparseMatrix_alloc: %ld x %ld system matrix has been allocated.\n",(long)out->numRows,(long)out->numCols);
    #endif
    return out;
  }
}

/* returns a reference to Paso_SparseMatrix in */

Paso_SparseMatrix* Paso_SparseMatrix_getReference(Paso_SparseMatrix* in) {
   if (in!=NULL) ++(in->reference_counter);
   return in;
}

/* deallocates a SparseMatrix: */

void Paso_SparseMatrix_free(Paso_SparseMatrix* in) {
  if (in!=NULL) {
     in->reference_counter--;
     if (in->reference_counter<=0) {
        Paso_Pattern_free(in->pattern);
        MEMFREE(in->val);
        MEMFREE(in);
        #ifdef Paso_TRACE
        printf("Paso_SparseMatrix_free: system matrix as been deallocated.\n");
        #endif
     }
   }
}

Paso_SparseMatrix* Paso_SparseMatrix_loadMM_toCSR( char *fileName_p )
{
    index_t *col_ind = NULL;
    index_t *row_ind = NULL;
    index_t *row_ptr = NULL;
    double *val = NULL;
    FILE *fileHandle_p = NULL;
        Paso_Pattern* mainPattern=NULL;
    Paso_SparseMatrix *out = NULL;
    int i, curr_row, scan_ret;
    MM_typecode matrixCode;
        Paso_resetError();

    /* open the file */
    fileHandle_p = fopen( fileName_p, "r" );
    if( fileHandle_p == NULL )
    {
        Paso_setError(IO_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Cannot read file for reading.");
        return NULL;
    }

    /* process banner */
    if( mm_read_banner(fileHandle_p, &matrixCode) != 0 )
    {
        Paso_setError(IO_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Error processing MM banner.");
        fclose( fileHandle_p );
        return NULL;
    }
    if( !(mm_is_real(matrixCode) && mm_is_sparse(matrixCode) && mm_is_general(matrixCode)) )
    {

        Paso_setError(TYPE_ERROR,"Paso_SparseMatrix_loadMM_toCSR: found Matrix Market type is not supported.");
        fclose( fileHandle_p );
        return NULL;
    }

    /* get matrix size */
    if( mm_read_mtx_crd_size(fileHandle_p, &M, &N, &nz) != 0 )
    {
        Paso_setError(IO_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Could not parse matrix size");
        fclose( fileHandle_p );
        return NULL;
    }

    /* prepare storage */
    col_ind = MEMALLOC( nz, index_t );
    row_ind = MEMALLOC( nz, index_t );
    val = MEMALLOC( nz, double );

    row_ptr = MEMALLOC( (M+1), index_t );

    if( col_ind == NULL || row_ind == NULL || val == NULL || row_ptr == NULL )
    {
        Paso_setError(MEMORY_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Could not allocate memory" );
        fclose( fileHandle_p );
        return NULL;
    }

    /* perform actual read of elements */
    for( i=0; i<nz; i++ )
    {
        scan_ret = fscanf( fileHandle_p, "%d %d %le\n", &row_ind[i], &col_ind[i], &val[i] );
        if (scan_ret!=3)
        {
            MEMFREE( val );
            MEMFREE( row_ind );
            MEMFREE( col_ind );
            MEMFREE( row_ptr );
            fclose(fileHandle_p);
            return NULL;
        }
        row_ind[i]--;
        col_ind[i]--;
    }
    fclose( fileHandle_p );

    /* sort the entries */
    q_sort( row_ind, col_ind, val, 0, nz );

    /* setup row_ptr */
    curr_row = 0;
    for( i=0; (i<nz && curr_row<M); curr_row++ )
    {
        while( row_ind[i] != curr_row )
            i++;
        row_ptr[curr_row] = i;
    }
    row_ptr[M] = nz;

        mainPattern=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,1,1,M,N,row_ptr,col_ind);
    out  = Paso_SparseMatrix_alloc(MATRIX_FORMAT_DEFAULT, mainPattern, 1, 1);
    /* copy values and cleanup temps */
    for( i=0; i<nz; i++ ) out->val[i] = val[i];
    
    Paso_Pattern_free(mainPattern);
    MEMFREE( val );
    MEMFREE( row_ind );
    return out;
}


void Paso_SparseMatrix_saveMM(Paso_SparseMatrix * A_p, char * fileName_p) {
  FILE * fileHandle_p = NULL;
  dim_t N,M,i, iptr_ij;
  MM_typecode matcode;                        

  if (A_p->col_block_size !=A_p->row_block_size) {
    Paso_setError(TYPE_ERROR, "Paso_SparseMatrix_saveMM: currently only square blocks are supported.");
    return;
  }
  if (A_p->row_block_size>3) {
       Paso_setError(TYPE_ERROR,"Paso_SparseMatrix_saveMM: currently only block size 3 is supported.\n");
       return;
  }

  if (A_p->type & MATRIX_FORMAT_SYM) {
    Paso_setError(TYPE_ERROR,"Paso_SparseMatrix_saveMM does not support symmetric storage scheme");
    return;
  }
  
  /* open the file */
  fileHandle_p = fopen(fileName_p, "w");
  if (fileHandle_p==NULL) {
    Paso_setError(IO_ERROR,"file could not be opened for writing");
    return;
  }

  if (A_p->type & MATRIX_FORMAT_CSC) {
    Paso_setError(TYPE_ERROR,"Paso_SparseMatrix_saveMM does not support CSC yet.");
  } else {
    mm_initialize_typecode(&matcode);
    mm_set_matrix(&matcode);
    mm_set_coordinate(&matcode);
    mm_set_real(&matcode);

    N= A_p->numRows; 
    M= A_p->numCols; 
    mm_write_banner(fileHandle_p, matcode); 
    mm_write_mtx_crd_size(fileHandle_p, N, M, A_p->pattern->ptr[N]);
    if(A_p->row_block_size==1)
      for (i=0; i<N; i++) 
        for (iptr_ij=A_p->pattern->ptr[i];iptr_ij<A_p->pattern->ptr[i+1]; ++iptr_ij) {
          fprintf(fileHandle_p, "%d %d %25.15e\n", i+1, A_p->pattern->index[iptr_ij]+1, A_p->val[iptr_ij]);
        }

    if(A_p->row_block_size==2)
      for (i=0; i<N; i++) {
        for (iptr_ij=A_p->pattern->ptr[i];iptr_ij<A_p->pattern->ptr[i+1]; ++iptr_ij) {
          fprintf(fileHandle_p, "%d %d %25.15e\n", i*2+1, 2*(A_p->pattern->index[iptr_ij])+1, A_p->val[iptr_ij*4]);
          fprintf(fileHandle_p, "%d %d %25.15e\n", i*2+1, 2*(A_p->pattern->index[iptr_ij])+1+1, A_p->val[iptr_ij*4+2]);
          fprintf(fileHandle_p, "%d %d %25.15e\n", i*2+1+1, 2*(A_p->pattern->index[iptr_ij])+1, A_p->val[iptr_ij*4+1]);
          fprintf(fileHandle_p, "%d %d %25.15e\n", i*2+1+1, 2*(A_p->pattern->index[iptr_ij])+1+1, A_p->val[iptr_ij*4+3]);
        }
      }

    if(A_p->row_block_size==3)
      for (i=0; i<N; i++) {
       for (iptr_ij=A_p->pattern->ptr[i];iptr_ij<A_p->pattern->ptr[i+1]; ++iptr_ij) {
          fprintf(fileHandle_p, "%d %d %25.15e\n", i*3+1, 3*(A_p->pattern->index[iptr_ij])+1, A_p->val[iptr_ij*9]);
          fprintf(fileHandle_p, "%d %d %25.15e\n", i*3+1, 3*(A_p->pattern->index[iptr_ij])+1+1, A_p->val[iptr_ij*9+3]);
          fprintf(fileHandle_p, "%d %d %25.15e\n", i*3+1, 3*(A_p->pattern->index[iptr_ij])+2+1, A_p->val[iptr_ij*9+6]);
          fprintf(fileHandle_p, "%d %d %25.15e\n", i*3+1+1, 3*(A_p->pattern->index[iptr_ij])+1, A_p->val[iptr_ij*9+1]);
          fprintf(fileHandle_p, "%d %d %25.15e\n", i*3+1+1, 3*(A_p->pattern->index[iptr_ij])+1+1, A_p->val[iptr_ij*9+4]);
          fprintf(fileHandle_p, "%d %d %25.15e\n", i*3+1+1, 3*(A_p->pattern->index[iptr_ij])+2+1, A_p->val[iptr_ij*9+7]);
          fprintf(fileHandle_p, "%d %d %25.15e\n", i*3+2+1, 3*(A_p->pattern->index[iptr_ij])+1, A_p->val[iptr_ij*9+2]);
          fprintf(fileHandle_p, "%d %d %25.15e\n", i*3+2+1, 3*(A_p->pattern->index[iptr_ij])+1+1, A_p->val[iptr_ij*9+5]);
          fprintf(fileHandle_p, "%d %d %25.15e\n", i*3+2+1, 3*(A_p->pattern->index[iptr_ij])+2+1, A_p->val[iptr_ij*9+8]);
        }
      } 
     
  }

  /* close the file */
  fclose(fileHandle_p);
  
  return;
}

1	ksteube	1313
2			/*******************************************************
3	ksteube	1811	*
4			* Copyright (c) 2003-2008 by University of Queensland
5			* Earth Systems Science Computational Center (ESSCC)
6			* http://www.uq.edu.au/esscc
7			*
8			* Primary Business: Queensland, Australia
9			* Licensed under the Open Software License version 3.0
10			* http://www.opensource.org/licenses/osl-3.0.php
11			*
12			*******************************************************/
13	ksteube	1313
14	ksteube	1811
15	ksteube	1313	/**************************************************************/
16
17			/* Paso: SparseMatrix */
18
19			/**************************************************************/
20
21			/* Author: gross@access.edu.au */
22
23			/**************************************************************/
24
25			#include "Paso.h"
26			#include "SparseMatrix.h"
27			#include "TRILINOS.h"
28	artak	2240	#include "mmio.h"
29	ksteube	1313
30			/**************************************************************/
31	artak	2240	static void swap( index_t, index_t, double*, int, int );
32			static void q_sort( index_t, index_t, double*, int, int );
33			/static void print_entries( index_t, index_t, double );*/
34	ksteube	1313
35	artak	2240	static int M, N, nz;
36
37
38			/* debug: print the entries */
39			/*
40			void print_entries( index_t r, index_t c, double *v )
41			{
42			int i;
43
44			for( i=0; i<nz; i++ )
45			{
46			printf( "(%ld, %ld) == %e\n", (long)r[i], (long)c[i], v[i] );
47			}
48			}
49			*/
50
51			/* swap function */
52			void swap( index_t r, index_t c, double *v, int left, int right )
53			{
54			double v_temp;
55			index_t temp;
56
57			temp = r[left];
58			r[left] = r[right];
59			r[right] = temp;
60
61			temp = c[left];
62			c[left] = c[right];
63			c[right] = temp;
64
65			v_temp = v[left];
66			v[left] = v[right];
67			v[right] = v_temp;
68			}
69
70			void q_sort( index_t row, index_t col, double *val, int begin, int end )
71			{
72			int l, r;
73			index_t pivot, lval;
74
75			if( end > begin )
76			{
77			pivot = N * row[begin] + col[begin];
78			l = begin + 1;
79			r = end;
80
81			while( l < r )
82			{
83			lval = N * row[l] + col[l];
84			if( lval < pivot )
85			l++;
86			else
87			{
88			r--;
89			swap( row, col, val, l, r );
90			}
91			}
92			l--;
93			swap( row, col, val, begin, l );
94			q_sort( row, col, val, begin, l );
95			q_sort( row, col, val, r, end );
96			}
97			}
98
99
100	ksteube	1313	/* allocates a SparseMatrix of type type using the given matrix pattern
101			if type is UNKOWN CSR is used.
102			if CSC or CSC_BLK1 is used pattern has to give the CSC pattern.
103			if CSR or CSR_BLK1 is used pattern has to give the CSR pattern.
104			Values are initialized by zero. */
105			Paso_SparseMatrix* Paso_SparseMatrix_alloc(Paso_SparseMatrixType type,Paso_Pattern *pattern, int row_block_size, int col_block_size) {
106
107			double time0;
108			Paso_SparseMatrix*out=NULL;
109			Paso_SparseMatrixType pattern_format_out;
110
111			Paso_resetError();
112			time0=Paso_timer();
113			out=MEMALLOC(1,Paso_SparseMatrix);
114			if (! Paso_checkPtr(out)) {
115			out->pattern=NULL;
116			out->val=NULL;
117			out->reference_counter=1;
118			out->type=type;
119
120			pattern_format_out= (type & MATRIX_FORMAT_OFFSET1)? PATTERN_FORMAT_OFFSET1: PATTERN_FORMAT_DEFAULT;
121			/* ====== compressed sparse columns === */
122			if (type & MATRIX_FORMAT_CSC) {
123			if (type & MATRIX_FORMAT_SYM) {
124			Paso_setError(TYPE_ERROR,"Generation of matrix pattern for symmetric CSC is not implemented yet.");
125			return NULL;
126			} else {
127			if ((type & MATRIX_FORMAT_BLK1) \|\| row_block_size!=col_block_size \|\| col_block_size>3) {
128			out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,col_block_size,row_block_size);
129			out->row_block_size=1;
130			out->col_block_size=1;
131			} else {
132			out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,1,1);
133			out->row_block_size=row_block_size;
134			out->col_block_size=col_block_size;
135			}
136			out->numRows = out->pattern->numInput;
137			out->numCols = out->pattern->numOutput;
138			}
139			} else {
140			/* ====== compressed sparse row === */
141			if (type & MATRIX_FORMAT_SYM) {
142			Paso_setError(TYPE_ERROR,"Generation of matrix pattern for symmetric CSR is not implemented yet.");
143			return NULL;
144			} else {
145			if ((type & MATRIX_FORMAT_BLK1) \|\| row_block_size!=col_block_size \|\| col_block_size>3) {
146			out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,row_block_size,col_block_size);
147			out->row_block_size=1;
148			out->col_block_size=1;
149			} else {
150			out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,1,1);
151			out->row_block_size=row_block_size;
152			out->col_block_size=col_block_size;
153			}
154			out->numRows = out->pattern->numOutput;
155			out->numCols = out->pattern->numInput;
156			}
157			}
158			out->block_size=out->row_block_size*out->col_block_size;
159			out->len=(size_t)(out->pattern->len)*(size_t)(out->block_size);
160
161			out->val=MEMALLOC(out->len,double);
162			if (! Paso_checkPtr(out->val)) Paso_SparseMatrix_setValues(out,DBLE(0));
163			}
164			/* all done: */
165			if (! Paso_noError()) {
166			Paso_SparseMatrix_free(out);
167			return NULL;
168			} else {
169			#ifdef Paso_TRACE
170			printf("timing: system matrix %.4e sec\n",Paso_timer()-time0);
171			printf("Paso_SparseMatrix_alloc: %ld x %ld system matrix has been allocated.\n",(long)out->numRows,(long)out->numCols);
172			#endif
173			return out;
174			}
175			}
176
177			/* returns a reference to Paso_SparseMatrix in */
178
179			Paso_SparseMatrix* Paso_SparseMatrix_getReference(Paso_SparseMatrix* in) {
180			if (in!=NULL) ++(in->reference_counter);
181			return in;
182			}
183
184			/* deallocates a SparseMatrix: */
185
186			void Paso_SparseMatrix_free(Paso_SparseMatrix* in) {
187			if (in!=NULL) {
188			in->reference_counter--;
189			if (in->reference_counter<=0) {
190			Paso_Pattern_free(in->pattern);
191			MEMFREE(in->val);
192			MEMFREE(in);
193			#ifdef Paso_TRACE
194			printf("Paso_SparseMatrix_free: system matrix as been deallocated.\n");
195			#endif
196			}
197			}
198			}
199	artak	2240
200			Paso_SparseMatrix* Paso_SparseMatrix_loadMM_toCSR( char *fileName_p )
201			{
202			index_t *col_ind = NULL;
203			index_t *row_ind = NULL;
204			index_t *row_ptr = NULL;
205			double *val = NULL;
206			FILE *fileHandle_p = NULL;
207			Paso_Pattern* mainPattern=NULL;
208			Paso_SparseMatrix *out = NULL;
209			int i, curr_row, scan_ret;
210			MM_typecode matrixCode;
211			Paso_resetError();
212
213			/* open the file */
214			fileHandle_p = fopen( fileName_p, "r" );
215			if( fileHandle_p == NULL )
216			{
217			Paso_setError(IO_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Cannot read file for reading.");
218			return NULL;
219			}
220
221			/* process banner */
222			if( mm_read_banner(fileHandle_p, &matrixCode) != 0 )
223			{
224			Paso_setError(IO_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Error processing MM banner.");
225			fclose( fileHandle_p );
226			return NULL;
227			}
228			if( !(mm_is_real(matrixCode) && mm_is_sparse(matrixCode) && mm_is_general(matrixCode)) )
229			{
230
231			Paso_setError(TYPE_ERROR,"Paso_SparseMatrix_loadMM_toCSR: found Matrix Market type is not supported.");
232			fclose( fileHandle_p );
233			return NULL;
234			}
235
236			/* get matrix size */
237			if( mm_read_mtx_crd_size(fileHandle_p, &M, &N, &nz) != 0 )
238			{
239			Paso_setError(IO_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Could not parse matrix size");
240			fclose( fileHandle_p );
241			return NULL;
242			}
243
244			/* prepare storage */
245			col_ind = MEMALLOC( nz, index_t );
246			row_ind = MEMALLOC( nz, index_t );
247			val = MEMALLOC( nz, double );
248
249			row_ptr = MEMALLOC( (M+1), index_t );
250
251			if( col_ind == NULL \|\| row_ind == NULL \|\| val == NULL \|\| row_ptr == NULL )
252			{
253			Paso_setError(MEMORY_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Could not allocate memory" );
254			fclose( fileHandle_p );
255			return NULL;
256			}
257
258			/* perform actual read of elements */
259			for( i=0; i<nz; i++ )
260			{
261			scan_ret = fscanf( fileHandle_p, "%d %d %le\n", &row_ind[i], &col_ind[i], &val[i] );
262			if (scan_ret!=3)
263			{
264			MEMFREE( val );
265			MEMFREE( row_ind );
266			MEMFREE( col_ind );
267			MEMFREE( row_ptr );
268			fclose(fileHandle_p);
269			return NULL;
270			}
271			row_ind[i]--;
272			col_ind[i]--;
273			}
274			fclose( fileHandle_p );
275
276			/* sort the entries */
277			q_sort( row_ind, col_ind, val, 0, nz );
278
279			/* setup row_ptr */
280			curr_row = 0;
281			for( i=0; (i<nz && curr_row<M); curr_row++ )
282			{
283			while( row_ind[i] != curr_row )
284			i++;
285			row_ptr[curr_row] = i;
286			}
287			row_ptr[M] = nz;
288
289			mainPattern=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,1,1,M,N,row_ptr,col_ind);
290			out = Paso_SparseMatrix_alloc(MATRIX_FORMAT_DEFAULT, mainPattern, 1, 1);
291			/* copy values and cleanup temps */
292			for( i=0; i<nz; i++ ) out->val[i] = val[i];
293
294			Paso_Pattern_free(mainPattern);
295			MEMFREE( val );
296			MEMFREE( row_ind );
297			return out;
298			}
299
300
301			void Paso_SparseMatrix_saveMM(Paso_SparseMatrix * A_p, char * fileName_p) {
302			FILE * fileHandle_p = NULL;
303			dim_t N,M,i, iptr_ij;
304			MM_typecode matcode;
305
306			if (A_p->col_block_size !=A_p->row_block_size) {
307			Paso_setError(TYPE_ERROR, "Paso_SparseMatrix_saveMM: currently only square blocks are supported.");
308			return;
309			}
310			if (A_p->row_block_size>3) {
311			Paso_setError(TYPE_ERROR,"Paso_SparseMatrix_saveMM: currently only block size 3 is supported.\n");
312			return;
313			}
314
315			if (A_p->type & MATRIX_FORMAT_SYM) {
316			Paso_setError(TYPE_ERROR,"Paso_SparseMatrix_saveMM does not support symmetric storage scheme");
317			return;
318			}
319
320			/* open the file */
321			fileHandle_p = fopen(fileName_p, "w");
322			if (fileHandle_p==NULL) {
323			Paso_setError(IO_ERROR,"file could not be opened for writing");
324			return;
325			}
326
327			if (A_p->type & MATRIX_FORMAT_CSC) {
328			Paso_setError(TYPE_ERROR,"Paso_SparseMatrix_saveMM does not support CSC yet.");
329			} else {
330			mm_initialize_typecode(&matcode);
331			mm_set_matrix(&matcode);
332			mm_set_coordinate(&matcode);
333			mm_set_real(&matcode);
334
335			N= A_p->numRows;
336			M= A_p->numCols;
337			mm_write_banner(fileHandle_p, matcode);
338			mm_write_mtx_crd_size(fileHandle_p, N, M, A_p->pattern->ptr[N]);
339			if(A_p->row_block_size==1)
340			for (i=0; i<N; i++)
341			for (iptr_ij=A_p->pattern->ptr[i];iptr_ij<A_p->pattern->ptr[i+1]; ++iptr_ij) {
342			fprintf(fileHandle_p, "%d %d %25.15e\n", i+1, A_p->pattern->index[iptr_ij]+1, A_p->val[iptr_ij]);
343			}
344
345			if(A_p->row_block_size==2)
346			for (i=0; i<N; i++) {
347			for (iptr_ij=A_p->pattern->ptr[i];iptr_ij<A_p->pattern->ptr[i+1]; ++iptr_ij) {
348			fprintf(fileHandle_p, "%d %d %25.15e\n", i2+1, 2(A_p->pattern->index[iptr_ij])+1, A_p->val[iptr_ij*4]);
349			fprintf(fileHandle_p, "%d %d %25.15e\n", i2+1, 2(A_p->pattern->index[iptr_ij])+1+1, A_p->val[iptr_ij*4+2]);
350			fprintf(fileHandle_p, "%d %d %25.15e\n", i2+1+1, 2(A_p->pattern->index[iptr_ij])+1, A_p->val[iptr_ij*4+1]);
351			fprintf(fileHandle_p, "%d %d %25.15e\n", i2+1+1, 2(A_p->pattern->index[iptr_ij])+1+1, A_p->val[iptr_ij*4+3]);
352			}
353			}
354
355			if(A_p->row_block_size==3)
356			for (i=0; i<N; i++) {
357			for (iptr_ij=A_p->pattern->ptr[i];iptr_ij<A_p->pattern->ptr[i+1]; ++iptr_ij) {
358			fprintf(fileHandle_p, "%d %d %25.15e\n", i3+1, 3(A_p->pattern->index[iptr_ij])+1, A_p->val[iptr_ij*9]);
359			fprintf(fileHandle_p, "%d %d %25.15e\n", i3+1, 3(A_p->pattern->index[iptr_ij])+1+1, A_p->val[iptr_ij*9+3]);
360			fprintf(fileHandle_p, "%d %d %25.15e\n", i3+1, 3(A_p->pattern->index[iptr_ij])+2+1, A_p->val[iptr_ij*9+6]);
361			fprintf(fileHandle_p, "%d %d %25.15e\n", i3+1+1, 3(A_p->pattern->index[iptr_ij])+1, A_p->val[iptr_ij*9+1]);
362			fprintf(fileHandle_p, "%d %d %25.15e\n", i3+1+1, 3(A_p->pattern->index[iptr_ij])+1+1, A_p->val[iptr_ij*9+4]);
363			fprintf(fileHandle_p, "%d %d %25.15e\n", i3+1+1, 3(A_p->pattern->index[iptr_ij])+2+1, A_p->val[iptr_ij*9+7]);
364			fprintf(fileHandle_p, "%d %d %25.15e\n", i3+2+1, 3(A_p->pattern->index[iptr_ij])+1, A_p->val[iptr_ij*9+2]);
365			fprintf(fileHandle_p, "%d %d %25.15e\n", i3+2+1, 3(A_p->pattern->index[iptr_ij])+1+1, A_p->val[iptr_ij*9+5]);
366			fprintf(fileHandle_p, "%d %d %25.15e\n", i3+2+1, 3(A_p->pattern->index[iptr_ij])+2+1, A_p->val[iptr_ij*9+8]);
367			}
368			}
369
370			}
371
372			/* close the file */
373			fclose(fileHandle_p);
374
375			return;
376			}
377