1 #include "adamat.h" /*I "mat.h" I*/ 2 3 #undef __FUNCT__ 4 #define __FUNCT__ "MatCreateADA" 5 /*@C 6 MatCreateADA - Creates a matrix M=A^T D1 A + D2 where D1, D2 are diagonal 7 8 Collective on matrix 9 10 Input Parameters: 11 + mat - matrix of arbitrary type 12 . d1 - A vector with diagonal elements of D1 13 - d2 - A vector 14 15 Output Parameters: 16 . J - New matrix whose operations are defined in terms of mat, D1, and D2. 17 18 Notes: 19 The user provides the input data and is responsible for destroying 20 this data after matrix J has been destroyed. 21 The operation MatMult(A,D2,D1) must be well defined. 22 Before calling the operation MatGetDiagonal(), the function 23 MatADAComputeDiagonal() must be called. The matrices A and D1 must 24 be the same during calls to MatADAComputeDiagonal() and 25 MatGetDiagonal(). 26 27 Level: developer 28 29 .seealso: MatCreate() 30 @*/ 31 PetscErrorCode MatCreateADA(Mat mat,Vec d1, Vec d2, Mat *J) 32 { 33 MPI_Comm comm=((PetscObject)mat)->comm; 34 TaoMatADACtx ctx; 35 PetscErrorCode ierr; 36 PetscInt nloc,n; 37 38 PetscFunctionBegin; 39 40 ierr = PetscNew(&ctx);CHKERRQ(ierr); 41 42 ctx->A=mat; 43 ctx->D1=d1; 44 ctx->D2=d2; 45 if (d1){ 46 ierr = VecDuplicate(d1,&ctx->W);CHKERRQ(ierr); 47 ierr = PetscObjectReference((PetscObject)d1);CHKERRQ(ierr); 48 } else { 49 ctx->W=0; 50 } 51 if (d2){ 52 ierr = VecDuplicate(d2,&ctx->W2);CHKERRQ(ierr); 53 ierr = VecDuplicate(d2,&ctx->ADADiag);CHKERRQ(ierr); 54 ierr = PetscObjectReference((PetscObject)d2);CHKERRQ(ierr); 55 } else { 56 ctx->W2=0; 57 ctx->ADADiag=0; 58 } 59 60 ctx->GotDiag=0; 61 ierr = PetscObjectReference((PetscObject)mat);CHKERRQ(ierr); 62 63 ierr=VecGetLocalSize(d2,&nloc);CHKERRQ(ierr); 64 ierr=VecGetSize(d2,&n);CHKERRQ(ierr); 65 66 ierr = MatCreateShell(comm,nloc,nloc,n,n,ctx,J);CHKERRQ(ierr); 67 68 ierr = MatShellSetOperation(*J,MATOP_MULT,(void(*)(void))MatMult_ADA);CHKERRQ(ierr); 69 ierr = MatShellSetOperation(*J,MATOP_DESTROY,(void(*)(void))MatDestroy_ADA);CHKERRQ(ierr); 70 ierr = MatShellSetOperation(*J,MATOP_VIEW,(void(*)(void))MatView_ADA);CHKERRQ(ierr); 71 ierr = MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_ADA);CHKERRQ(ierr); 72 ierr = MatShellSetOperation(*J,MATOP_DIAGONAL_SET,(void(*)(void))MatDiagonalSet_ADA);CHKERRQ(ierr); 73 ierr = MatShellSetOperation(*J,MATOP_SHIFT,(void(*)(void))MatShift_ADA);CHKERRQ(ierr); 74 ierr = MatShellSetOperation(*J,MATOP_EQUAL,(void(*)(void))MatEqual_ADA);CHKERRQ(ierr); 75 ierr = MatShellSetOperation(*J,MATOP_SCALE,(void(*)(void))MatScale_ADA);CHKERRQ(ierr); 76 ierr = MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)(void))MatTranspose_ADA);CHKERRQ(ierr); 77 ierr = MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_ADA);CHKERRQ(ierr); 78 ierr = MatShellSetOperation(*J,MATOP_GET_SUBMATRICES,(void(*)(void))MatGetSubMatrices_ADA);CHKERRQ(ierr); 79 ierr = MatShellSetOperation(*J,MATOP_NORM,(void(*)(void))MatNorm_ADA);CHKERRQ(ierr); 80 ierr = MatShellSetOperation(*J,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_ADA);CHKERRQ(ierr); 81 ierr = MatShellSetOperation(*J,MATOP_GET_SUBMATRIX,(void(*)(void))MatGetSubMatrix_ADA);CHKERRQ(ierr); 82 83 ierr = PetscLogObjectParent((PetscObject)(*J),(PetscObject)ctx->W);CHKERRQ(ierr); 84 ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)(*J));CHKERRQ(ierr); 85 86 ierr = MatSetOption(*J,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 87 PetscFunctionReturn(0); 88 } 89 90 #undef __FUNCT__ 91 #define __FUNCT__ "MatMult_ADA" 92 PetscErrorCode MatMult_ADA(Mat mat,Vec a,Vec y) 93 { 94 TaoMatADACtx ctx; 95 PetscReal one = 1.0; 96 PetscErrorCode ierr; 97 98 PetscFunctionBegin; 99 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 100 101 ierr = MatMult(ctx->A,a,ctx->W);CHKERRQ(ierr); 102 if (ctx->D1){ 103 ierr = VecPointwiseMult(ctx->W,ctx->D1,ctx->W);CHKERRQ(ierr); 104 } 105 ierr = MatMultTranspose(ctx->A,ctx->W,y);CHKERRQ(ierr); 106 if (ctx->D2){ 107 ierr = VecPointwiseMult(ctx->W2, ctx->D2, a);CHKERRQ(ierr); 108 ierr = VecAXPY(y, one, ctx->W2);CHKERRQ(ierr); 109 } 110 PetscFunctionReturn(0); 111 } 112 113 #undef __FUNCT__ 114 #define __FUNCT__ "MatMultTranspose_ADA" 115 PetscErrorCode MatMultTranspose_ADA(Mat mat,Vec a,Vec y) 116 { 117 PetscErrorCode ierr; 118 119 PetscFunctionBegin; 120 ierr = MatMult_ADA(mat,a,y);CHKERRQ(ierr); 121 PetscFunctionReturn(0); 122 } 123 124 #undef __FUNCT__ 125 #define __FUNCT__ "MatDiagonalSet_ADA" 126 PetscErrorCode MatDiagonalSet_ADA(Vec D, Mat M) 127 { 128 TaoMatADACtx ctx; 129 PetscReal zero=0.0,one = 1.0; 130 PetscErrorCode ierr; 131 132 PetscFunctionBegin; 133 ierr = MatShellGetContext(M,(void **)&ctx);CHKERRQ(ierr); 134 135 if (ctx->D2==NULL){ 136 ierr = VecDuplicate(D,&ctx->D2);CHKERRQ(ierr); 137 ierr = VecSet(ctx->D2, zero);CHKERRQ(ierr); 138 } 139 ierr = VecAXPY(ctx->D2, one, D);CHKERRQ(ierr); 140 141 PetscFunctionReturn(0); 142 } 143 144 #undef __FUNCT__ 145 #define __FUNCT__ "MatDestroy_ADA" 146 PetscErrorCode MatDestroy_ADA(Mat mat) 147 { 148 PetscErrorCode ierr; 149 TaoMatADACtx ctx; 150 151 PetscFunctionBegin; 152 ierr=MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 153 if (ctx->W) { 154 ierr=VecDestroy(&ctx->W);CHKERRQ(ierr); 155 } 156 if (ctx->W2) { 157 ierr=VecDestroy(&ctx->W2);CHKERRQ(ierr); 158 } 159 if (ctx->ADADiag) { 160 ierr=VecDestroy(&ctx->ADADiag);CHKERRQ(ierr); 161 } 162 if (ctx->A) { 163 ierr=MatDestroy(&ctx->A);CHKERRQ(ierr); 164 } 165 if (ctx->D1) { 166 ierr=VecDestroy(&ctx->D1);CHKERRQ(ierr); 167 } 168 if (ctx->D2) { 169 ierr=VecDestroy(&ctx->D2);CHKERRQ(ierr); 170 } 171 ierr = PetscFree(ctx);CHKERRQ(ierr); 172 PetscFunctionReturn(0); 173 } 174 175 #undef __FUNCT__ 176 #define __FUNCT__ "MatView_ADA" 177 PetscErrorCode MatView_ADA(Mat mat,PetscViewer viewer) 178 { 179 180 PetscFunctionBegin; 181 PetscFunctionReturn(0); 182 } 183 184 #undef __FUNCT__ 185 #define __FUNCT__ "MatShift_ADA" 186 PetscErrorCode MatShift_ADA(Mat Y, PetscReal a) 187 { 188 PetscErrorCode ierr; 189 TaoMatADACtx ctx; 190 191 PetscFunctionBegin; 192 ierr = MatShellGetContext(Y,(void **)&ctx);CHKERRQ(ierr); 193 ierr = VecShift(ctx->D2,a);CHKERRQ(ierr); 194 PetscFunctionReturn(0); 195 } 196 197 #undef __FUNCT__ 198 #define __FUNCT__ "MatDuplicate_ADA" 199 PetscErrorCode MatDuplicate_ADA(Mat mat,MatDuplicateOption op,Mat *M) 200 { 201 PetscErrorCode ierr; 202 TaoMatADACtx ctx; 203 Mat A2; 204 Vec D1b=NULL,D2b; 205 206 PetscFunctionBegin; 207 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 208 ierr = MatDuplicate(ctx->A,op,&A2);CHKERRQ(ierr); 209 if (ctx->D1){ 210 ierr = VecDuplicate(ctx->D1,&D1b);CHKERRQ(ierr); 211 ierr = VecCopy(ctx->D1,D1b);CHKERRQ(ierr); 212 } 213 ierr = VecDuplicate(ctx->D2,&D2b);CHKERRQ(ierr); 214 ierr = VecCopy(ctx->D2,D2b);CHKERRQ(ierr); 215 ierr = MatCreateADA(A2,D1b,D2b,M);CHKERRQ(ierr); 216 if (ctx->D1){ 217 ierr = PetscObjectDereference((PetscObject)D1b);CHKERRQ(ierr); 218 } 219 ierr = PetscObjectDereference((PetscObject)D2b);CHKERRQ(ierr); 220 ierr = PetscObjectDereference((PetscObject)A2);CHKERRQ(ierr); 221 222 PetscFunctionReturn(0); 223 } 224 225 #undef __FUNCT__ 226 #define __FUNCT__ "MatEqual_ADA" 227 PetscErrorCode MatEqual_ADA(Mat A,Mat B,PetscBool *flg) 228 { 229 PetscErrorCode ierr; 230 TaoMatADACtx ctx1,ctx2; 231 232 PetscFunctionBegin; 233 ierr = MatShellGetContext(A,(void **)&ctx1);CHKERRQ(ierr); 234 ierr = MatShellGetContext(B,(void **)&ctx2);CHKERRQ(ierr); 235 ierr = VecEqual(ctx1->D2,ctx2->D2,flg);CHKERRQ(ierr); 236 if (*flg==PETSC_TRUE){ 237 ierr = VecEqual(ctx1->D1,ctx2->D1,flg);CHKERRQ(ierr); 238 } 239 if (*flg==PETSC_TRUE){ 240 ierr = MatEqual(ctx1->A,ctx2->A,flg);CHKERRQ(ierr); 241 } 242 PetscFunctionReturn(0); 243 } 244 245 #undef __FUNCT__ 246 #define __FUNCT__ "MatScale_ADA" 247 PetscErrorCode MatScale_ADA(Mat mat, PetscReal a) 248 { 249 PetscErrorCode ierr; 250 TaoMatADACtx ctx; 251 252 PetscFunctionBegin; 253 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 254 ierr = VecScale(ctx->D1,a);CHKERRQ(ierr); 255 if (ctx->D2){ 256 ierr = VecScale(ctx->D2,a);CHKERRQ(ierr); 257 } 258 PetscFunctionReturn(0); 259 } 260 261 #undef __FUNCT__ 262 #define __FUNCT__ "MatTranspose_ADA" 263 PetscErrorCode MatTranspose_ADA(Mat mat,Mat *B) 264 { 265 PetscErrorCode ierr; 266 TaoMatADACtx ctx; 267 268 PetscFunctionBegin; 269 if (*B){ 270 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 271 ierr = MatDuplicate(mat,MAT_COPY_VALUES,B);CHKERRQ(ierr); 272 } 273 PetscFunctionReturn(0); 274 } 275 276 #undef __FUNCT__ 277 #define __FUNCT__ "MatADAComputeDiagonal" 278 PetscErrorCode MatADAComputeDiagonal(Mat mat) 279 { 280 PetscErrorCode ierr; 281 PetscInt i,m,n,low,high; 282 PetscReal *dtemp,*dptr; 283 TaoMatADACtx ctx; 284 285 PetscFunctionBegin; 286 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 287 288 ierr = MatGetOwnershipRange(mat, &low, &high);CHKERRQ(ierr); 289 ierr = MatGetSize(mat,&m,&n);CHKERRQ(ierr); 290 291 ierr = PetscMalloc( n*sizeof(PetscReal),&dtemp );CHKERRQ(ierr); 292 293 for (i=0; i<n; i++){ 294 ierr = MatGetColumnVector(ctx->A, ctx->W, i);CHKERRQ(ierr); 295 ierr = VecPointwiseMult(ctx->W,ctx->W,ctx->W);CHKERRQ(ierr); 296 ierr = VecDotBegin(ctx->D1, ctx->W,dtemp+i);CHKERRQ(ierr); 297 } 298 for (i=0; i<n; i++){ 299 ierr = VecDotEnd(ctx->D1, ctx->W,dtemp+i);CHKERRQ(ierr); 300 } 301 302 ierr = VecGetArray(ctx->ADADiag,&dptr);CHKERRQ(ierr); 303 for (i=low; i<high; i++){ 304 dptr[i-low]= dtemp[i]; 305 } 306 ierr = VecRestoreArray(ctx->ADADiag,&dptr);CHKERRQ(ierr); 307 if (dtemp) { 308 ierr = PetscFree(dtemp);CHKERRQ(ierr); 309 } 310 PetscFunctionReturn(0); 311 } 312 313 #undef __FUNCT__ 314 #define __FUNCT__ "MatGetDiagonal_ADA" 315 PetscErrorCode MatGetDiagonal_ADA(Mat mat,Vec v) 316 { 317 PetscErrorCode ierr; 318 PetscReal one=1.0; 319 TaoMatADACtx ctx; 320 321 PetscFunctionBegin; 322 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 323 ierr = MatADAComputeDiagonal(mat);CHKERRQ(ierr); 324 ierr=VecCopy(ctx->ADADiag,v);CHKERRQ(ierr); 325 if (ctx->D2){ 326 ierr=VecAXPY(v, one, ctx->D2);CHKERRQ(ierr); 327 } 328 329 PetscFunctionReturn(0); 330 } 331 332 #undef __FUNCT__ 333 #define __FUNCT__ "MatGetSubMatrices_ADA" 334 PetscErrorCode MatGetSubMatrices_ADA(Mat A,PetscInt n, IS *irow,IS *icol,MatReuse scall,Mat **B) 335 { 336 PetscErrorCode ierr; 337 PetscInt i; 338 339 PetscFunctionBegin; 340 if (scall == MAT_INITIAL_MATRIX) { 341 ierr = PetscMalloc( (n+1)*sizeof(Mat),B );CHKERRQ(ierr); 342 } 343 344 for ( i=0; i<n; i++ ) { 345 ierr = MatGetSubMatrix_ADA(A,irow[i],icol[i],scall,&(*B)[i]);CHKERRQ(ierr); 346 } 347 PetscFunctionReturn(0); 348 } 349 350 #undef __FUNCT__ 351 #define __FUNCT__ "MatGetSubMatrix_ADA" 352 PetscErrorCode MatGetSubMatrix_ADA(Mat mat,IS isrow,IS iscol,MatReuse cll, 353 Mat *newmat) 354 { 355 PetscErrorCode ierr; 356 PetscInt low,high; 357 PetscInt n,nlocal,i; 358 const PetscInt *iptr; 359 PetscReal *dptr,*ddptr,zero=0.0; 360 VecType type_name; 361 IS ISrow; 362 Vec D1,D2; 363 Mat Atemp; 364 TaoMatADACtx ctx; 365 366 PetscFunctionBegin; 367 368 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 369 370 ierr = MatGetOwnershipRange(ctx->A,&low,&high);CHKERRQ(ierr); 371 ierr = ISCreateStride(((PetscObject)mat)->comm,high-low,low,1,&ISrow);CHKERRQ(ierr); 372 ierr = MatGetSubMatrix(ctx->A,ISrow,iscol,cll,&Atemp);CHKERRQ(ierr); 373 ierr = ISDestroy(&ISrow);CHKERRQ(ierr); 374 375 if (ctx->D1){ 376 ierr=VecDuplicate(ctx->D1,&D1);CHKERRQ(ierr); 377 ierr=VecCopy(ctx->D1,D1);CHKERRQ(ierr); 378 } else { 379 D1=NULL; 380 } 381 382 if (ctx->D2){ 383 ierr=ISGetSize(isrow,&n);CHKERRQ(ierr); 384 ierr=ISGetLocalSize(isrow,&nlocal);CHKERRQ(ierr); 385 ierr=VecCreate(((PetscObject)(ctx->D2))->comm,&D2);CHKERRQ(ierr); 386 ierr=VecGetType(ctx->D2,&type_name);CHKERRQ(ierr); 387 ierr=VecSetSizes(D2,nlocal,n);CHKERRQ(ierr); 388 ierr=VecSetType(D2,type_name);CHKERRQ(ierr); 389 ierr=VecSet(D2, zero);CHKERRQ(ierr); 390 ierr=VecGetArray(ctx->D2, &dptr);CHKERRQ(ierr); 391 ierr=VecGetArray(D2, &ddptr);CHKERRQ(ierr); 392 ierr=ISGetIndices(isrow,&iptr);CHKERRQ(ierr); 393 for (i=0;i<nlocal;i++){ 394 ddptr[i] = dptr[iptr[i]-low]; 395 } 396 ierr=ISRestoreIndices(isrow,&iptr);CHKERRQ(ierr); 397 ierr=VecRestoreArray(D2, &ddptr);CHKERRQ(ierr); 398 ierr=VecRestoreArray(ctx->D2, &dptr);CHKERRQ(ierr); 399 400 } else { 401 D2=NULL; 402 } 403 404 ierr = MatCreateADA(Atemp,D1,D2,newmat);CHKERRQ(ierr); 405 ierr = MatShellGetContext(*newmat,(void **)&ctx);CHKERRQ(ierr); 406 ierr = PetscObjectDereference((PetscObject)Atemp);CHKERRQ(ierr); 407 if (ctx->D1){ 408 ierr = PetscObjectDereference((PetscObject)D1);CHKERRQ(ierr); 409 } 410 if (ctx->D2){ 411 ierr = PetscObjectDereference((PetscObject)D2);CHKERRQ(ierr); 412 } 413 PetscFunctionReturn(0); 414 } 415 416 #undef __FUNCT__ 417 #define __FUNCT__ "MatGetRowADA" 418 PetscErrorCode MatGetRowADA(Mat mat,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscReal **vals) 419 { 420 PetscErrorCode ierr; 421 PetscInt m,n; 422 423 PetscFunctionBegin; 424 ierr = MatGetSize(mat,&m,&n);CHKERRQ(ierr); 425 426 if (*ncols>0){ 427 ierr = PetscMalloc( (*ncols)*sizeof(PetscInt),cols );CHKERRQ(ierr); 428 ierr = PetscMalloc( (*ncols)*sizeof(PetscReal),vals );CHKERRQ(ierr); 429 } else { 430 *cols=NULL; 431 *vals=NULL; 432 } 433 434 PetscFunctionReturn(0); 435 } 436 437 #undef __FUNCT__ 438 #define __FUNCT__ "MatRestoreRowADA" 439 PetscErrorCode MatRestoreRowADA(Mat mat,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscReal **vals) 440 { 441 PetscErrorCode ierr; 442 PetscFunctionBegin; 443 if (*ncols>0){ 444 ierr = PetscFree(*cols); CHKERRQ(ierr); 445 ierr = PetscFree(*vals); CHKERRQ(ierr); 446 } 447 *cols=NULL; 448 *vals=NULL; 449 PetscFunctionReturn(0); 450 } 451 452 #undef __FUNCT__ 453 #define __FUNCT__ "MatGetColumnVectorADA" 454 PetscErrorCode MatGetColumnVector_ADA(Mat mat,Vec Y, PetscInt col) 455 { 456 PetscErrorCode ierr; 457 PetscInt low,high; 458 PetscReal zero=0.0,one=1.0; 459 460 PetscFunctionBegin; 461 ierr=VecSet(Y, zero);CHKERRQ(ierr); 462 ierr=VecGetOwnershipRange(Y,&low,&high);CHKERRQ(ierr); 463 if (col>=low && col<high){ 464 ierr=VecSetValue(Y,col,one,INSERT_VALUES);CHKERRQ(ierr); 465 } 466 ierr=VecAssemblyBegin(Y);CHKERRQ(ierr); 467 ierr=VecAssemblyEnd(Y);CHKERRQ(ierr); 468 ierr=MatMult_ADA(mat,Y,Y);CHKERRQ(ierr); 469 470 PetscFunctionReturn(0); 471 } 472 473 PetscErrorCode MatConvert_ADA(Mat mat,MatType newtype,Mat *NewMat) 474 { 475 PetscErrorCode ierr; 476 PetscMPIInt size; 477 PetscBool sametype, issame, isdense, isseqdense; 478 TaoMatADACtx ctx; 479 480 PetscFunctionBegin; 481 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 482 ierr = MPI_Comm_size(((PetscObject)mat)->comm,&size);CHKERRQ(ierr); 483 484 ierr = PetscObjectTypeCompare((PetscObject)mat,newtype,&sametype);CHKERRQ(ierr); 485 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSAME,&issame);CHKERRQ(ierr); 486 ierr = PetscObjectTypeCompare((PetscObject)mat,MATMPIDENSE,&isdense);CHKERRQ(ierr); 487 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSEQDENSE,&isseqdense);CHKERRQ(ierr); 488 489 if (sametype || issame) { 490 ierr=MatDuplicate(mat,MAT_COPY_VALUES,NewMat);CHKERRQ(ierr); 491 } else if (isdense) { 492 PetscInt i,j,low,high,m,n,M,N; 493 PetscReal *dptr; 494 Vec X; 495 496 ierr = VecDuplicate(ctx->D2,&X);CHKERRQ(ierr); 497 ierr=MatGetSize(mat,&M,&N);CHKERRQ(ierr); 498 ierr=MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 499 ierr = MatCreateDense(((PetscObject)mat)->comm,m,m,N,N,NULL,NewMat);CHKERRQ(ierr); 500 ierr = MatGetOwnershipRange(*NewMat,&low,&high);CHKERRQ(ierr); 501 for (i=0;i<M;i++){ 502 ierr = MatGetColumnVector_ADA(mat,X,i);CHKERRQ(ierr); 503 ierr = VecGetArray(X,&dptr);CHKERRQ(ierr); 504 for (j=0; j<high-low; j++){ 505 ierr = MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);CHKERRQ(ierr); 506 } 507 ierr=VecRestoreArray(X,&dptr);CHKERRQ(ierr); 508 } 509 ierr=MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 510 ierr=MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 511 ierr = VecDestroy(&X);CHKERRQ(ierr); 512 } else if (isseqdense && size==1){ 513 PetscInt i,j,low,high,m,n,M,N; 514 PetscReal *dptr; 515 Vec X; 516 517 ierr = VecDuplicate(ctx->D2,&X);CHKERRQ(ierr); 518 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 519 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 520 ierr = MatCreateSeqDense(((PetscObject)mat)->comm,N,N,NULL,NewMat);CHKERRQ(ierr); 521 ierr = MatGetOwnershipRange(*NewMat,&low,&high);CHKERRQ(ierr); 522 for (i=0;i<M;i++){ 523 ierr = MatGetColumnVector_ADA(mat,X,i);CHKERRQ(ierr); 524 ierr = VecGetArray(X,&dptr);CHKERRQ(ierr); 525 for (j=0; j<high-low; j++){ 526 ierr = MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);CHKERRQ(ierr); 527 } 528 ierr=VecRestoreArray(X,&dptr);CHKERRQ(ierr); 529 } 530 ierr=MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 531 ierr=MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 532 ierr=VecDestroy(&X);CHKERRQ(ierr); 533 } else SETERRQ(PETSC_COMM_SELF,1,"No support to convert objects to that type"); 534 PetscFunctionReturn(0); 535 } 536 537 #undef __FUNCT__ 538 #define __FUNCT__ "MatNorm_ADA" 539 PetscErrorCode MatNorm_ADA(Mat mat,NormType type,PetscReal *norm) 540 { 541 PetscErrorCode ierr; 542 TaoMatADACtx ctx; 543 544 PetscFunctionBegin; 545 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 546 if (type == NORM_FROBENIUS) { 547 *norm = 1.0; 548 } else if (type == NORM_1 || type == NORM_INFINITY) { 549 *norm = 1.0; 550 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 551 PetscFunctionReturn(0); 552 } 553