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
2	/*******************************************************
3	*
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
14
15	/**************************************************************/
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	#include "mmio.h"
29
30	/**************************************************************/
31	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
35	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	/* 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
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