paso/src/SparseMatrix_AMGcomponents.c


/*******************************************************
*
* Copyright (c) 2003-2009 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_AMGcomponents */

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

/* Author: Artak Amirbekyan, artak@uq.edu.au */

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

#include "Paso.h"
#include "SparseMatrix.h"
#include "PasoUtil.h"

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

    Methods nessecary for AMG preconditioner

*/

/* Prolongation matrix is constracted by W_FC as P<- [W_FC]
                                                     [I_CC]
*/

Paso_SparseMatrix* Paso_SparseMatrix_getProlongation(Paso_SparseMatrix* W, index_t* mis_marker){
 
  Paso_Pattern *outpattern=NULL;
  Paso_SparseMatrix *out=NULL;
  index_t iptr,wptr;
  
  Paso_IndexList* index_list=NULL;
  dim_t n=W->numRows+W->numCols;
  dim_t i,j,k=0;

  index_list=TMPMEMALLOC(n,Paso_IndexList);
   if (! Paso_checkPtr(index_list)) {
        #pragma omp parallel for private(i) schedule(static)
        for(i=0;i<n;++i) {
             index_list[i].extension=NULL;
             index_list[i].n=0;
        }
    }
  
    for (i=0;i<n;++i) {
      if (mis_marker[i]) {
          for (iptr=W->pattern->ptr[k];iptr<W->pattern->ptr[k+1]; ++iptr) {
          j=W->pattern->index[iptr];
          Paso_IndexList_insertIndex(&(index_list[i]),j);
        }
        k++;
      }
      else {
          Paso_IndexList_insertIndex(&(index_list[i]),i-k);
      }
    }
    
    outpattern=Paso_IndexList_createPattern(0, n,index_list,0,W->numCols,0);
    out=Paso_SparseMatrix_alloc(W->type,outpattern,1,1,TRUE);
    
    k=0;
    for (i=0;i<out->numRows;++i) {
      if (mis_marker[i]) {
        wptr=W->pattern->ptr[k];
        for (iptr=out->pattern->ptr[i];iptr<out->pattern->ptr[i+1]; ++iptr) {
            out->val[iptr]=W->val[wptr];
            wptr++;
           }
        k++;   
      }
      else {
           for (iptr=out->pattern->ptr[i];iptr<out->pattern->ptr[i+1]; ++iptr){
             out->val[iptr]=1;
           }
      }
    }
    
     /* clean up */
   if (index_list!=NULL) {
        #pragma omp parallel for private(i) schedule(static)
        for(i=0;i<n;++i) Paso_IndexList_free(index_list[i].extension);
     }
    TMPMEMFREE(index_list);
    Paso_Pattern_free(outpattern);
    return out;
}


/* Restriction matrix R=P^T */

Paso_SparseMatrix* Paso_SparseMatrix_getRestriction(Paso_SparseMatrix* P){
 
  Paso_Pattern *outpattern=NULL;
  Paso_SparseMatrix *out=NULL;
  
  Paso_IndexList* index_list=NULL;
  dim_t C=P->numCols;
  dim_t F=P->numRows-C;
  dim_t n=C+F;
  dim_t i,j,k=0;
  index_t iptr,jptr;

  index_list=TMPMEMALLOC(C,Paso_IndexList);
   if (! Paso_checkPtr(index_list)) {
        #pragma omp parallel for private(i) schedule(static)
        for(i=0;i<C;++i) {
             index_list[i].extension=NULL;
             index_list[i].n=0;
        }
    }
  

    /*#pragma omp parallel for private(i,iptr,j) schedule(static)*/
    for (i=0;i<n;++i) {
          for (iptr=P->pattern->ptr[i];iptr<P->pattern->ptr[i+1]; ++iptr) {
             j=P->pattern->index[iptr];
             Paso_IndexList_insertIndex(&(index_list[j]),i);
        }
    }
   
    outpattern=Paso_IndexList_createPattern(0, C,index_list,0,C+F,0);
    out=Paso_SparseMatrix_alloc(P->type,outpattern,1,1,TRUE);
    
    for (i=0;i<out->numRows;++i) {
           for (iptr=out->pattern->ptr[i];iptr<out->pattern->ptr[i+1]; ++iptr) {
                 j=out->pattern->index[iptr];
                  /*This for later will be replaced by bsearch!!*/
                  for (jptr=P->pattern->ptr[j];jptr<P->pattern->ptr[j+1]; ++jptr) {
                        k=P->pattern->index[jptr];
                        if(k==i) {
                             out->val[iptr]=P->val[jptr];
                        }
                  }
           }
      }
    
     /* clean up */
   if (index_list!=NULL) {
        #pragma omp parallel for private(i) schedule(static)
        for(i=0;i<C;++i) Paso_IndexList_free(index_list[i].extension);
     }
    TMPMEMFREE(index_list);
    Paso_Pattern_free(outpattern);
    return out;
}


void Paso_SparseMatrix_updateWeights(Paso_SparseMatrix* A,Paso_SparseMatrix* W_FC, index_t* mis_marker){
 
  double *alpha;
  double *beta;
  double *alpha_den;
  double *beta_den;
  double a_ii=0;
  
  dim_t n=A->numRows;
  dim_t F=W_FC->numRows;
  dim_t i,j,k;
  
  index_t iPtr,*index, *where_p;
 
  alpha=TMPMEMALLOC(F,double);
  beta=TMPMEMALLOC(F,double);
  alpha_den=TMPMEMALLOC(F,double);
  beta_den=TMPMEMALLOC(F,double);

  k=0;  
  for (i = 0; i < n; ++i) {
      if(mis_marker[i]) {
            iPtr=A->pattern->ptr[i];
            index=&(A->pattern->index[iPtr]);
            where_p=(index_t*)bsearch(&i,
                                     index,
                                     A->pattern->ptr[i + 1]-A->pattern->ptr[i],
                                     sizeof(index_t),
                                     Paso_comparIndex);
            if (where_p!=NULL) {
                   a_ii=A->val[iPtr+(index_t)(where_p-index)];
            }
            else {
             Paso_setError(VALUE_ERROR,"Paso_Solver_getAMG: Main diagonal element of system matrix is missing.");
            }
            alpha[k]=0;
            beta[k]=0;
            alpha_den[k]=0;
            beta_den[k]=0;
            /*compute beta_i=sum_{j\inN_i}a^{+}_ij and alpha_i=sum_{j\inN_i}a^{-}_ij and denominators for both alpha_i and beta_i*/
            for (iPtr=A->pattern->ptr[i];iPtr<A->pattern->ptr[i + 1]; ++iPtr) {
                 j=A->pattern->index[iPtr];
                 if(j!=i) {
                        if(A->val[iPtr]<0) {
                          alpha[k]+=A->val[iPtr];
                          if(!mis_marker[j]) {
                            alpha_den[k]+=A->val[iPtr];
                          }
                        }
                        else if(A->val[iPtr]>0) {
                          beta[k]+=A->val[iPtr];
                          if(!mis_marker[j]) {
                            beta_den[k]+=A->val[iPtr];
                          }
                        }
                 }
             }
            if(alpha_den[k]!=0) {
                 alpha[k]=alpha[k]/(alpha_den[k]*a_ii);
            }
            if(beta_den[k]!=0) {
                 beta[k]=beta[k]/(beta_den[k]*a_ii);
            }
            k++;
       /*printf("Got in row=%d, alpha[%d]=%e, beta[%d]=%e, a_den=%e, b_den=%e \n",i,k-1,alpha[k-1],k-1,beta[k-1],alpha_den[k-1],beta_den[k-1]);*/
      }
   }
  
      for (i = 0; i < W_FC->numRows; ++i) {
            for (iPtr=W_FC->pattern->ptr[i];iPtr<W_FC->pattern->ptr[i + 1]; ++iPtr) {
                 j=W_FC->pattern->index[iPtr];
                 /*if(!mis_marker[j]) {*/
                   if(W_FC->val[iPtr]<0) {
                      W_FC->val[iPtr]=-alpha[i]*W_FC->val[iPtr];
                    }
                   else if(W_FC->val[iPtr]>0) {
                      W_FC->val[iPtr]=-beta[i]*W_FC->val[iPtr];
                   }
                 /*}*/
            }
         }
      
      TMPMEMFREE(alpha);
      TMPMEMFREE(beta);
      TMPMEMFREE(alpha_den);
      TMPMEMFREE(beta_den);
}

/*A_c=R*A*P taking into account coarseneing */
/*
void Paso_Solver_getCoarseMatrix(Paso_SparseMatrix* A_c, Paso_SparseMatrix* A,Paso_SparseMatrix *R,Paso_SparseMatrix *P) {

  index_t iptrA_c,iptrR,iptrP,iptrA;
  dim_t i,j,k,l,m;
  double first_sum=0,second_sum=0,p_lj=0,a_kl=0,r_ik=0;
 */
  /*a^c_ij=sum_k^n(r_ik)sum_l^n(a_kl*P_lj)*/
 /* 
  for(i = 0; i < A_c->numRows; i++) {
     for(iptrA_c = A_c->pattern->ptr[i]; iptrA_c < A_c->pattern->ptr[i+1]; ++iptrA_c) {
      j = A_c->pattern->index[iptrA_c];
      second_sum=0;
      for(iptrR = R->pattern->ptr[i]; iptrR < R->pattern->ptr[i+1]; ++iptrR) {
        k = R->pattern->index[iptrR];
        first_sum=0;
        for(iptrA = A->pattern->ptr[k]; iptrA < A->pattern->ptr[k+1]; ++iptrA) {
         l= A->pattern->index[iptrA];
 */        /*This loop has to be replaced by bsearch */
  /*       for(iptrP = P->pattern->ptr[l]; iptrP < P->pattern->ptr[l+1]; ++iptrP) {
           m=P->pattern->index[iptrP];
           if(m==j) {
            p_lj=P->val[iptrP];
            break;
           }
         }
         a_kl=A->val[iptrA];
         first_sum+=a_kl*p_lj;
        }
        r_ik=R->val[iptrR];
        second_sum+=r_ik*first_sum;
        if(i==0) {
        printf("row %d, col %d, k=%d, firstsum = %e, r_ik=%e \n",i,j,k,first_sum,r_ik);
        }
      }
      if(i==0) {
        printf("row %d, col %d, secondsum = %e \n",i,j,second_sum);
      }
      A_c->val[iptrA_c]=second_sum;
     }
    }

}
  */
  
/*A_c=R*A*P taking into account coarseneing */
/*Paso_SparseMatrix* Paso_Solver_getCoarseMatrix(Paso_SparseMatrix* A, Paso_SparseMatrix *R, Paso_SparseMatrix *P) {

Paso_SparseMatrix *temp=NULL;
Paso_SparseMatrix *out=NULL;

    temp=Paso_SparseMatrix_MatrixMatrix(A,P);
    out=Paso_SparseMatrix_MatrixMatrix(R,temp);
    
    Paso_SparseMatrix_free(temp);
      
return out;      
}
 */

Paso_SparseMatrix* Paso_SparseMatrix_MatrixMatrix(Paso_SparseMatrix* A, Paso_SparseMatrix* B) {
  index_t iptrA,iptrB,iptrC;
  dim_t i,j=0,k,l;
  Paso_Pattern* outpattern=NULL;
  Paso_SparseMatrix *out=NULL;
  double sum,b_lj=0;
  
  outpattern=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,A->pattern,B->pattern);
  out=Paso_SparseMatrix_alloc(A->type,outpattern,1,1, TRUE);
  
  
  for(i = 0; i < out->numRows; i++) {
    for(iptrC = out->pattern->ptr[i]; iptrC < out->pattern->ptr[i+1]; ++iptrC) {
      j = out->pattern->index[iptrC];
      sum=0;
      b_lj=0;
      for(iptrA = A->pattern->ptr[i]; iptrA < A->pattern->ptr[i+1]; ++iptrA) {
            k=A->pattern->index[iptrA];
            /*can be replaced by bsearch */
            for(iptrB = B->pattern->ptr[k]; iptrB < B->pattern->ptr[k+1]; ++iptrB) {
               l=B->pattern->index[iptrB];
               if(l==j) {
                  b_lj=B->val[iptrB];
                  break;
               }
            }
            sum+=A->val[iptrA]*b_lj;
      }
      out->val[iptrC]=sum;
   }
  }
      
  Paso_Pattern_free(outpattern);

return out;      
}


Paso_SparseMatrix* Paso_Solver_getCoarseMatrix(Paso_SparseMatrix* A,Paso_SparseMatrix *R,Paso_SparseMatrix *P) {

  index_t iptrA_c,iptrR,iptrP,iptrA;
  dim_t i,j,k,l,m;
  double first_sum=0,second_sum=0,p_lj=0,a_kl=0,r_ik=0;
  Paso_Pattern* temp=NULL;
  Paso_Pattern* outpattern=NULL;
  Paso_SparseMatrix *A_c=NULL;
  
  temp=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,A->pattern,P->pattern);
  outpattern=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,R->pattern,temp);
  A_c=Paso_SparseMatrix_alloc(A->type,outpattern,1,1, TRUE);
  
  
  /*a^c_ij=sum_k^n(r_ik)sum_l^n(a_kl*P_lj)*/

  for(i = 0; i < A_c->numRows; i++) {
     for(iptrA_c = A_c->pattern->ptr[i]; iptrA_c < A_c->pattern->ptr[i+1]; ++iptrA_c) {
      j = A_c->pattern->index[iptrA_c];
      second_sum=0;
      for(iptrR = R->pattern->ptr[i]; iptrR < R->pattern->ptr[i+1]; ++iptrR) {
        k = R->pattern->index[iptrR];
        first_sum=0;
        for(iptrA = A->pattern->ptr[k]; iptrA < A->pattern->ptr[k+1]; ++iptrA) {
         l= A->pattern->index[iptrA];
         p_lj=0;
        /*This loop has to be replaced by bsearch */
         for(iptrP = P->pattern->ptr[l]; iptrP < P->pattern->ptr[l+1]; ++iptrP) {
           m=P->pattern->index[iptrP];
           if(m==j) {
            p_lj=P->val[iptrP];
            break;
           }
         }
         a_kl=A->val[iptrA];
         first_sum+=a_kl*p_lj;
        }
        r_ik=R->val[iptrR];
        second_sum+=r_ik*first_sum;
      }
      A_c->val[iptrA_c]=second_sum;
     }
    }
    
    Paso_Pattern_free(outpattern);
    Paso_Pattern_free(temp);
    
return A_c;
}

1
2	/*******************************************************
3	*
4	* Copyright (c) 2003-2009 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_AMGcomponents */
18
19	/**************************************************************/
20
21	/* Author: Artak Amirbekyan, artak@uq.edu.au */
22
23	/**************************************************************/
24
25	#include "Paso.h"
26	#include "SparseMatrix.h"
27	#include "PasoUtil.h"
28
29	/**************************************************************
30
31	Methods nessecary for AMG preconditioner
32
33	*/
34
35	/* Prolongation matrix is constracted by W_FC as P<- [W_FC]
36	[I_CC]
37	*/
38
39	Paso_SparseMatrix* Paso_SparseMatrix_getProlongation(Paso_SparseMatrix* W, index_t* mis_marker){
40
41	Paso_Pattern *outpattern=NULL;
42	Paso_SparseMatrix *out=NULL;
43	index_t iptr,wptr;
44
45	Paso_IndexList* index_list=NULL;
46	dim_t n=W->numRows+W->numCols;
47	dim_t i,j,k=0;
48
49	index_list=TMPMEMALLOC(n,Paso_IndexList);
50	if (! Paso_checkPtr(index_list)) {
51	#pragma omp parallel for private(i) schedule(static)
52	for(i=0;i<n;++i) {
53	index_list[i].extension=NULL;
54	index_list[i].n=0;
55	}
56	}
57
58	for (i=0;i<n;++i) {
59	if (mis_marker[i]) {
60	for (iptr=W->pattern->ptr[k];iptr<W->pattern->ptr[k+1]; ++iptr) {
61	j=W->pattern->index[iptr];
62	Paso_IndexList_insertIndex(&(index_list[i]),j);
63	}
64	k++;
65	}
66	else {
67	Paso_IndexList_insertIndex(&(index_list[i]),i-k);
68	}
69	}
70
71	outpattern=Paso_IndexList_createPattern(0, n,index_list,0,W->numCols,0);
72	out=Paso_SparseMatrix_alloc(W->type,outpattern,1,1,TRUE);
73
74	k=0;
75	for (i=0;i<out->numRows;++i) {
76	if (mis_marker[i]) {
77	wptr=W->pattern->ptr[k];
78	for (iptr=out->pattern->ptr[i];iptr<out->pattern->ptr[i+1]; ++iptr) {
79	out->val[iptr]=W->val[wptr];
80	wptr++;
81	}
82	k++;
83	}
84	else {
85	for (iptr=out->pattern->ptr[i];iptr<out->pattern->ptr[i+1]; ++iptr){
86	out->val[iptr]=1;
87	}
88	}
89	}
90
91	/* clean up */
92	if (index_list!=NULL) {
93	#pragma omp parallel for private(i) schedule(static)
94	for(i=0;i<n;++i) Paso_IndexList_free(index_list[i].extension);
95	}
96	TMPMEMFREE(index_list);
97	Paso_Pattern_free(outpattern);
98	return out;
99	}
100
101
102	/* Restriction matrix R=P^T */
103
104	Paso_SparseMatrix* Paso_SparseMatrix_getRestriction(Paso_SparseMatrix* P){
105
106	Paso_Pattern *outpattern=NULL;
107	Paso_SparseMatrix *out=NULL;
108
109	Paso_IndexList* index_list=NULL;
110	dim_t C=P->numCols;
111	dim_t F=P->numRows-C;
112	dim_t n=C+F;
113	dim_t i,j,k=0;
114	index_t iptr,jptr;
115
116	index_list=TMPMEMALLOC(C,Paso_IndexList);
117	if (! Paso_checkPtr(index_list)) {
118	#pragma omp parallel for private(i) schedule(static)
119	for(i=0;i<C;++i) {
120	index_list[i].extension=NULL;
121	index_list[i].n=0;
122	}
123	}
124
125
126	/#pragma omp parallel for private(i,iptr,j) schedule(static)/
127	for (i=0;i<n;++i) {
128	for (iptr=P->pattern->ptr[i];iptr<P->pattern->ptr[i+1]; ++iptr) {
129	j=P->pattern->index[iptr];
130	Paso_IndexList_insertIndex(&(index_list[j]),i);
131	}
132	}
133
134	outpattern=Paso_IndexList_createPattern(0, C,index_list,0,C+F,0);
135	out=Paso_SparseMatrix_alloc(P->type,outpattern,1,1,TRUE);
136
137	for (i=0;i<out->numRows;++i) {
138	for (iptr=out->pattern->ptr[i];iptr<out->pattern->ptr[i+1]; ++iptr) {
139	j=out->pattern->index[iptr];
140	/This for later will be replaced by bsearch!!/
141	for (jptr=P->pattern->ptr[j];jptr<P->pattern->ptr[j+1]; ++jptr) {
142	k=P->pattern->index[jptr];
143	if(k==i) {
144	out->val[iptr]=P->val[jptr];
145	}
146	}
147	}
148	}
149
150	/* clean up */
151	if (index_list!=NULL) {
152	#pragma omp parallel for private(i) schedule(static)
153	for(i=0;i<C;++i) Paso_IndexList_free(index_list[i].extension);
154	}
155	TMPMEMFREE(index_list);
156	Paso_Pattern_free(outpattern);
157	return out;
158	}
159
160
161	void Paso_SparseMatrix_updateWeights(Paso_SparseMatrix* A,Paso_SparseMatrix* W_FC, index_t* mis_marker){
162
163	double *alpha;
164	double *beta;
165	double *alpha_den;
166	double *beta_den;
167	double a_ii=0;
168
169	dim_t n=A->numRows;
170	dim_t F=W_FC->numRows;
171	dim_t i,j,k;
172
173	index_t iPtr,index, where_p;
174
175	alpha=TMPMEMALLOC(F,double);
176	beta=TMPMEMALLOC(F,double);
177	alpha_den=TMPMEMALLOC(F,double);
178	beta_den=TMPMEMALLOC(F,double);
179
180	k=0;
181	for (i = 0; i < n; ++i) {
182	if(mis_marker[i]) {
183	iPtr=A->pattern->ptr[i];
184	index=&(A->pattern->index[iPtr]);
185	where_p=(index_t*)bsearch(&i,
186	index,
187	A->pattern->ptr[i + 1]-A->pattern->ptr[i],
188	sizeof(index_t),
189	Paso_comparIndex);
190	if (where_p!=NULL) {
191	a_ii=A->val[iPtr+(index_t)(where_p-index)];
192	}
193	else {
194	Paso_setError(VALUE_ERROR,"Paso_Solver_getAMG: Main diagonal element of system matrix is missing.");
195	}
196	alpha[k]=0;
197	beta[k]=0;
198	alpha_den[k]=0;
199	beta_den[k]=0;
200	/compute beta_i=sum_{j\inN_i}a^{+}_ij and alpha_i=sum_{j\inN_i}a^{-}_ij and denominators for both alpha_i and beta_i/
201	for (iPtr=A->pattern->ptr[i];iPtr<A->pattern->ptr[i + 1]; ++iPtr) {
202	j=A->pattern->index[iPtr];
203	if(j!=i) {
204	if(A->val[iPtr]<0) {
205	alpha[k]+=A->val[iPtr];
206	if(!mis_marker[j]) {
207	alpha_den[k]+=A->val[iPtr];
208	}
209	}
210	else if(A->val[iPtr]>0) {
211	beta[k]+=A->val[iPtr];
212	if(!mis_marker[j]) {
213	beta_den[k]+=A->val[iPtr];
214	}
215	}
216	}
217	}
218	if(alpha_den[k]!=0) {
219	alpha[k]=alpha[k]/(alpha_den[k]*a_ii);
220	}
221	if(beta_den[k]!=0) {
222	beta[k]=beta[k]/(beta_den[k]*a_ii);
223	}
224	k++;
225	/printf("Got in row=%d, alpha[%d]=%e, beta[%d]=%e, a_den=%e, b_den=%e \n",i,k-1,alpha[k-1],k-1,beta[k-1],alpha_den[k-1],beta_den[k-1]);/
226	}
227	}
228
229	for (i = 0; i < W_FC->numRows; ++i) {
230	for (iPtr=W_FC->pattern->ptr[i];iPtr<W_FC->pattern->ptr[i + 1]; ++iPtr) {
231	j=W_FC->pattern->index[iPtr];
232	/if(!mis_marker[j]) {/
233	if(W_FC->val[iPtr]<0) {
234	W_FC->val[iPtr]=-alpha[i]*W_FC->val[iPtr];
235	}
236	else if(W_FC->val[iPtr]>0) {
237	W_FC->val[iPtr]=-beta[i]*W_FC->val[iPtr];
238	}
239	/}/
240	}
241	}
242
243	TMPMEMFREE(alpha);
244	TMPMEMFREE(beta);
245	TMPMEMFREE(alpha_den);
246	TMPMEMFREE(beta_den);
247	}
248
249	/A_c=RAP taking into account coarseneing /
250	/*
251	void Paso_Solver_getCoarseMatrix(Paso_SparseMatrix* A_c, Paso_SparseMatrix* A,Paso_SparseMatrix R,Paso_SparseMatrix P) {
252
253	index_t iptrA_c,iptrR,iptrP,iptrA;
254	dim_t i,j,k,l,m;
255	double first_sum=0,second_sum=0,p_lj=0,a_kl=0,r_ik=0;
256	*/
257	/a^c_ij=sum_k^n(r_ik)sum_l^n(a_klP_lj)*/
258	/*
259	for(i = 0; i < A_c->numRows; i++) {
260	for(iptrA_c = A_c->pattern->ptr[i]; iptrA_c < A_c->pattern->ptr[i+1]; ++iptrA_c) {
261	j = A_c->pattern->index[iptrA_c];
262	second_sum=0;
263	for(iptrR = R->pattern->ptr[i]; iptrR < R->pattern->ptr[i+1]; ++iptrR) {
264	k = R->pattern->index[iptrR];
265	first_sum=0;
266	for(iptrA = A->pattern->ptr[k]; iptrA < A->pattern->ptr[k+1]; ++iptrA) {
267	l= A->pattern->index[iptrA];
268	/ /This loop has to be replaced by bsearch */
269	/* for(iptrP = P->pattern->ptr[l]; iptrP < P->pattern->ptr[l+1]; ++iptrP) {
270	m=P->pattern->index[iptrP];
271	if(m==j) {
272	p_lj=P->val[iptrP];
273	break;
274	}
275	}
276	a_kl=A->val[iptrA];
277	first_sum+=a_kl*p_lj;
278	}
279	r_ik=R->val[iptrR];
280	second_sum+=r_ik*first_sum;
281	if(i==0) {
282	printf("row %d, col %d, k=%d, firstsum = %e, r_ik=%e \n",i,j,k,first_sum,r_ik);
283	}
284	}
285	if(i==0) {
286	printf("row %d, col %d, secondsum = %e \n",i,j,second_sum);
287	}
288	A_c->val[iptrA_c]=second_sum;
289	}
290	}
291
292	}
293	*/
294
295	/A_c=RAP taking into account coarseneing /
296	/Paso_SparseMatrix Paso_Solver_getCoarseMatrix(Paso_SparseMatrix* A, Paso_SparseMatrix R, Paso_SparseMatrix P) {
297
298	Paso_SparseMatrix *temp=NULL;
299	Paso_SparseMatrix *out=NULL;
300
301	temp=Paso_SparseMatrix_MatrixMatrix(A,P);
302	out=Paso_SparseMatrix_MatrixMatrix(R,temp);
303
304	Paso_SparseMatrix_free(temp);
305
306	return out;
307	}
308	*/
309
310	Paso_SparseMatrix* Paso_SparseMatrix_MatrixMatrix(Paso_SparseMatrix* A, Paso_SparseMatrix* B) {
311	index_t iptrA,iptrB,iptrC;
312	dim_t i,j=0,k,l;
313	Paso_Pattern* outpattern=NULL;
314	Paso_SparseMatrix *out=NULL;
315	double sum,b_lj=0;
316
317	outpattern=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,A->pattern,B->pattern);
318	out=Paso_SparseMatrix_alloc(A->type,outpattern,1,1, TRUE);
319
320
321	for(i = 0; i < out->numRows; i++) {
322	for(iptrC = out->pattern->ptr[i]; iptrC < out->pattern->ptr[i+1]; ++iptrC) {
323	j = out->pattern->index[iptrC];
324	sum=0;
325	b_lj=0;
326	for(iptrA = A->pattern->ptr[i]; iptrA < A->pattern->ptr[i+1]; ++iptrA) {
327	k=A->pattern->index[iptrA];
328	/can be replaced by bsearch /
329	for(iptrB = B->pattern->ptr[k]; iptrB < B->pattern->ptr[k+1]; ++iptrB) {
330	l=B->pattern->index[iptrB];
331	if(l==j) {
332	b_lj=B->val[iptrB];
333	break;
334	}
335	}
336	sum+=A->val[iptrA]*b_lj;
337	}
338	out->val[iptrC]=sum;
339	}
340	}
341
342	Paso_Pattern_free(outpattern);
343
344	return out;
345	}
346
347
348	Paso_SparseMatrix* Paso_Solver_getCoarseMatrix(Paso_SparseMatrix* A,Paso_SparseMatrix R,Paso_SparseMatrix P) {
349
350	index_t iptrA_c,iptrR,iptrP,iptrA;
351	dim_t i,j,k,l,m;
352	double first_sum=0,second_sum=0,p_lj=0,a_kl=0,r_ik=0;
353	Paso_Pattern* temp=NULL;
354	Paso_Pattern* outpattern=NULL;
355	Paso_SparseMatrix *A_c=NULL;
356
357	temp=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,A->pattern,P->pattern);
358	outpattern=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,R->pattern,temp);
359	A_c=Paso_SparseMatrix_alloc(A->type,outpattern,1,1, TRUE);
360
361
362	/a^c_ij=sum_k^n(r_ik)sum_l^n(a_klP_lj)*/
363
364	for(i = 0; i < A_c->numRows; i++) {
365	for(iptrA_c = A_c->pattern->ptr[i]; iptrA_c < A_c->pattern->ptr[i+1]; ++iptrA_c) {
366	j = A_c->pattern->index[iptrA_c];
367	second_sum=0;
368	for(iptrR = R->pattern->ptr[i]; iptrR < R->pattern->ptr[i+1]; ++iptrR) {
369	k = R->pattern->index[iptrR];
370	first_sum=0;
371	for(iptrA = A->pattern->ptr[k]; iptrA < A->pattern->ptr[k+1]; ++iptrA) {
372	l= A->pattern->index[iptrA];
373	p_lj=0;
374	/This loop has to be replaced by bsearch /
375	for(iptrP = P->pattern->ptr[l]; iptrP < P->pattern->ptr[l+1]; ++iptrP) {
376	m=P->pattern->index[iptrP];
377	if(m==j) {
378	p_lj=P->val[iptrP];
379	break;
380	}
381	}
382	a_kl=A->val[iptrA];
383	first_sum+=a_kl*p_lj;
384	}
385	r_ik=R->val[iptrR];
386	second_sum+=r_ik*first_sum;
387	}
388	A_c->val[iptrA_c]=second_sum;
389	}
390	}
391
392	Paso_Pattern_free(outpattern);
393	Paso_Pattern_free(temp);
394
395	return A_c;
396	}
397