source: issm/trunk/src/c/Reducematrixfromgtofx/Reducematrixfromgton.cpp@ 3223

/*!\file: Reducematrixfromgton
 * \brief
 */ 

#ifdef HAVE_CONFIG_H
        #include "config.h"
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#undef __FUNCT__ 
#define __FUNCT__ "Reducematrixfromgton"

#include "../toolkits/toolkits.h"


void  Reducematrixfromgton(Mat* pKnn,Mat Kgg,Mat Gmn,double* pv_m,int msize, double* pv_n,int nsize,int flag){


        /*output: */
        Mat Knn=NULL;

        /*intermediary: */
        Mat Kmm=NULL;
        Mat Kmn=NULL;
        Mat Knm=NULL;
        Mat KnmGmn=NULL;
        Mat KmmGmn=NULL;
        Mat tGmn=NULL;
        Mat tGmnKmn=NULL;
        PetscScalar a=1;

        if (msize){
                
                /*Partition K_gg*/
                MatPartition(&Kmm, Kgg, pv_m,msize,pv_m,msize); 
                MatPartition(&Kmn, Kgg, pv_m,msize,pv_n,nsize); 
                MatPartition(&Knm, Kgg, pv_n,nsize,pv_m,msize); 
                MatPartition(&Knn, Kgg, pv_n,nsize,pv_n,nsize);
        
                /*Reduce K_gg to K_nn*/
                //K_nn = K_nn + K_nm * G_mn;
                MatMatMult(Knm,Gmn,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&KnmGmn);
                MatAXPY(Knn,a,KnmGmn,DIFFERENT_NONZERO_PATTERN);

                if (flag!=1){
                        //K_mn = K_mn + K_mm * G_mn;
                        MatMatMult(Kmm,Gmn,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&KmmGmn);
                        MatAXPY(Kmn,a,KmmGmn,DIFFERENT_NONZERO_PATTERN);
                        
                        //K_nn = K_nn + G_mn' * K_mn;
                        #if _PETSC_VERSION_ == 2
                        MatTranspose(Gmn,&tGmn);
                        #else
                        MatTranspose(Gmn,MAT_INITIAL_MATRIX,&tGmn);
                        #endif
                        MatMatMult(tGmn,Kmn,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&tGmnKmn);
                        MatAXPY(Knn,a,tGmnKmn,DIFFERENT_NONZERO_PATTERN);
                }
        }
        else{
                /*msize=0 just return input*/
                MatDuplicate(Kgg,MAT_COPY_VALUES,&Knn);
        }       

        /*Free ressources:*/
        MatFree(&Kmm);
        MatFree(&Kmn);
        MatFree(&Knm);
        MatFree(&KnmGmn);
        MatFree(&KmmGmn);
        MatFree(&tGmn);
        MatFree(&tGmnKmn);
        
        /*Assign output pointers:*/
        *pKnn=Knn;
}