lmvm/dense/denseqn.c

640   Mat_DQN  *lbfgs = (Mat_DQN *)lmvm->ctx;
644   if (!lbfgs->temp_mat) PetscCall(MatDuplicate(lbfgs->YtS_triu_strict, MAT_SHARE_NONZERO_PATTERN, &lbfgs->temp_mat));
645   PetscCall(MatCopy(lbfgs->YtS_triu_strict, lbfgs->temp_mat, SAME_NONZERO_PATTERN));
646   PetscCall(MatDiagonalScale(lbfgs->temp_mat, lbfgs->inv_diag_vec, NULL));
649   PetscCall(MatConjugate(lbfgs->temp_mat));
652     PetscCall(MatDenseGetLocalMatrix(lbfgs->YtS_triu_strict, &YtS_local));
653     PetscCall(MatDenseGetLocalMatrix(lbfgs->temp_mat, &temp_local));
664   Mat_DQN  *lbfgs = (Mat_DQN *)lmvm->ctx;
675   if (!lbfgs->YtS_triu_strict) {
676     PetscCall(MatDuplicate(lbfgs->StY_triu, MAT_SHARE_NONZERO_PATTERN, &lbfgs->YtS_triu_strict));
677     PetscCall(MatDestroy(&lbfgs->StBS));
678     PetscCall(MatDuplicate(lbfgs->StY_triu, MAT_SHARE_NONZERO_PATTERN, &lbfgs->StBS));
679     PetscCall(MatDestroy(&lbfgs->J));
680     PetscCall(MatDuplicate(lbfgs->StY_triu, MAT_SHARE_NONZERO_PATTERN, &lbfgs->J));
681     PetscCall(MatDestroy(&lbfgs->BS));
682     PetscCall(MatDuplicate(Yfull, MAT_SHARE_NONZERO_PATTERN, &lbfgs->BS));
683     PetscCall(MatShift(lbfgs->StBS, 1.0));
684     lbfgs->num_mult_updates = oldest_update(m, k);
686   if (lbfgs->num_mult_updates == k) PetscFunctionReturn(PETSC_SUCCESS);
694     PetscInt StBS_idx = lbfgs->strategy == MAT_LMVM_DENSE_INPLACE ? S_idx : history_index(m, k, j);
696     PetscCall(MatDenseGetColumnVecWrite(lbfgs->BS, S_idx, &Bs_j));
700     PetscCall(MatDenseGetColumnVecWrite(lbfgs->StBS, StBS_idx, &StBs_j));
702     lbfgs->St_count++;
703     if (lbfgs->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecRecycleOrderToHistoryOrder(B, StBs_j, lbfgs->num_updates, lbfgs->cyclic_work_vec));
704     PetscCall(MatDenseRestoreColumnVecWrite(lbfgs->StBS, StBS_idx, &StBs_j));
705     PetscCall(MatDenseRestoreColumnVecWrite(lbfgs->BS, S_idx, &Bs_j));
707   prev_oldest = oldest_update(m, lbfgs->num_mult_updates);
708   if (lbfgs->strategy == MAT_LMVM_DENSE_REORDER && prev_oldest < oldest_update(m, k)) {
712     PetscCall(MatMove_LR3(B, lbfgs->YtS_triu_strict, m_keep));
714   PetscCall(MatGetLocalSize(lbfgs->YtS_triu_strict, &m_local, NULL));
715   j_0 = PetscMax(lbfgs->num_mult_updates, oldest_update(m, k));
718     PetscInt YtS_idx = lbfgs->strategy == MAT_LMVM_DENSE_INPLACE ? S_idx : history_index(m, k, j);
722     PetscCall(MatDenseGetColumnVecWrite(lbfgs->YtS_triu_strict, YtS_idx, &Yts_j));
724     lbfgs->Yt_count++;
725     if (lbfgs->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecRecycleOrderToHistoryOrder(B, Yts_j, lbfgs->num_updates, lbfgs->cyclic_work_vec));
726     PetscCall(MatDenseRestoreColumnVecWrite(lbfgs->YtS_triu_strict, YtS_idx, &Yts_j));
732       PetscCall(MatDenseGetLocalMatrix(lbfgs->YtS_triu_strict, &YtS_local));
738   if (!lbfgs->inv_diag_vec) PetscCall(VecDuplicate(lbfgs->diag_vec, &lbfgs->inv_diag_vec));
739   PetscCall(VecCopy(lbfgs->diag_vec, lbfgs->inv_diag_vec));
740   PetscCall(VecReciprocal(lbfgs->inv_diag_vec));
741   PetscCall(MatDenseGetLocalMatrix(lbfgs->J, &J_local));
743   PetscCall(MatGetLDLT(B, lbfgs->J));
744   PetscCall(MatAXPY(lbfgs->J, 1.0, lbfgs->StBS, SAME_NONZERO_PATTERN));
749   lbfgs->num_mult_updates = lbfgs->num_updates;
761   Mat_DQN         *lbfgs  = (Mat_DQN *)lmvm->ctx;
762   Vec              rwork1 = lbfgs->rwork1;
775   if (!lbfgs->num_updates) {
781   if (F == lbfgs->Fprev_ref && Fstate == lbfgs->Fprev_state) {
782     PetscCall(VecCopy(lbfgs->StFprev, rwork1));
785     lbfgs->St_count++;
789   if (lbfgs->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecRecycleOrderToHistoryOrder(H, rwork1, lbfgs->num_updates, lbfgs->cyclic_work_vec));
790   PetscCall(MatUpperTriangularSolveInPlace(H, lbfgs->StY_triu, rwork1, PETSC_FALSE, lbfgs->num_updates, lbfgs->strategy));
792   if (lbfgs->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecHistoryOrderToRecycleOrder(H, rwork1, lbfgs->num_updates, lbfgs->cyclic_work_vec));
794   PetscCall(VecCopy(F, lbfgs->column_work));
795   PetscCall(MatMultAddColumnRange(Yfull, rwork1, lbfgs->column_work, lbfgs->column_work, 0, h));
796   lbfgs->Y_count++;
798   PetscCall(VecPointwiseMult(rwork1, lbfgs->diag_vec_recycle_order, rwork1));
799   PetscCall(MatDQNApplyJ0Inv(H, lbfgs->column_work, dX));
802   lbfgs->Yt_count++;
804   if (lbfgs->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecRecycleOrderToHistoryOrder(H, rwork1, lbfgs->num_updates, lbfgs->cyclic_work_vec));
805   PetscCall(MatUpperTriangularSolveInPlace(H, lbfgs->StY_triu, rwork1, PETSC_TRUE, lbfgs->num_updates, lbfgs->strategy));
807   if (lbfgs->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecHistoryOrderToRecycleOrder(H, rwork1, lbfgs->num_updates, lbfgs->cyclic_work_vec));
810   lbfgs->S_count++;
849   Mat_DQN  *lbfgs = (Mat_DQN *)lmvm->ctx;
865   if (!lbfgs->num_updates) PetscFunctionReturn(PETSC_SUCCESS); /* No updates stored yet */
867   if (lbfgs->use_recursive) {
875     PetscCall(MatMultHermitianTransposeColumnRange(Yfull, X, lbfgs->rwork1, 0, h));
876     lbfgs->Yt_count++;
878     PetscCall(VecGetLocalSize(lbfgs->rwork1, &local_n));
880     if (lbfgs->needPQ) {
885         PetscCall(MatGetColumnVector(Sfull, lbfgs->column_work, idx));
886         PetscCall(MatDQNApplyJ0Fwd(B, lbfgs->column_work, lbfgs->PQ[idx]));
887         PetscCall(MatMultHermitianTransposeColumnRange(Yfull, lbfgs->column_work, lbfgs->rwork3, 0, h));
888         PetscCall(VecGetArrayAndMemType(lbfgs->rwork3, &workscalar, &memtype));
903           PetscCall(MatGetColumnVector(Sfull, lbfgs->column_work2, idx_j));
904           PetscCall(VecDot(lbfgs->PQ[idx], lbfgs->column_work2, &sjtpi));
906           PetscCall(MatGetColumnVector(Yfull, lbfgs->column_work2, idx_j));
908           PetscCall(VecAXPBYPCZ(lbfgs->PQ[idx], -sjtpi / lbfgs->stp[idx_j], yjtsi / lbfgs->yts[idx_j], 1.0, lbfgs->PQ[idx_j], lbfgs->column_work2));
910         PetscCall(VecDot(lbfgs->PQ[idx], lbfgs->column_work, &stp));
911         lbfgs->stp[idx] = PetscRealPart(stp);
913       lbfgs->needPQ = PETSC_FALSE;
916     PetscCall(VecGetArrayAndMemType(lbfgs->rwork1, &workscalar, &memtype));
931       PetscCall(MatGetColumnVector(Sfull, lbfgs->column_work, idx));
932       PetscCall(MatGetColumnVector(Yfull, lbfgs->column_work2, idx));
933       PetscCall(VecDot(Z, lbfgs->column_work, &stz));
934       PetscCall(VecAXPBYPCZ(Z, -stz / lbfgs->stp[idx], ytx / lbfgs->yts[idx], 1.0, lbfgs->PQ[idx], lbfgs->column_work2));
936     PetscCall(VecRestoreArrayAndMemType(lbfgs->rwork1, &workscalar));
939     PetscCall(MatMultHermitianTransposeColumnRange(Yfull, X, lbfgs->rwork1, 0, h));
940     lbfgs->Yt_count++;
941     PetscCall(MatMultHermitianTransposeColumnRange(Sfull, Z, lbfgs->rwork2, 0, h));
942     lbfgs->St_count++;
943     if (lbfgs->strategy == MAT_LMVM_DENSE_REORDER) {
944       PetscCall(VecRecycleOrderToHistoryOrder(B, lbfgs->rwork1, lbfgs->num_updates, lbfgs->cyclic_work_vec));
945       PetscCall(VecRecycleOrderToHistoryOrder(B, lbfgs->rwork2, lbfgs->num_updates, lbfgs->cyclic_work_vec));
948     PetscCall(VecPointwiseMult(lbfgs->rwork1, lbfgs->rwork1, lbfgs->inv_diag_vec));
949     if (PetscDefined(USE_COMPLEX)) PetscCall(MatConjugate(lbfgs->YtS_triu_strict));
950     PetscCall(MatMultTransposeAdd(lbfgs->YtS_triu_strict, lbfgs->rwork1, lbfgs->rwork2, lbfgs->rwork2));
951     if (PetscDefined(USE_COMPLEX)) PetscCall(MatConjugate(lbfgs->YtS_triu_strict));
953     if (!lbfgs->rwork2_local) PetscCall(VecCreateLocalVector(lbfgs->rwork2, &lbfgs->rwork2_local));
954     if (!lbfgs->rwork3_local) PetscCall(VecCreateLocalVector(lbfgs->rwork3, &lbfgs->rwork3_local));
955     PetscCall(VecGetLocalVectorRead(lbfgs->rwork2, lbfgs->rwork2_local));
956     PetscCall(VecGetLocalVector(lbfgs->rwork3, lbfgs->rwork3_local));
957     PetscCall(MatDenseGetLocalMatrix(lbfgs->J, &J_local));
958     PetscCall(VecGetSize(lbfgs->rwork2_local, &m_local));
960       PetscCall(MatDenseGetLocalMatrix(lbfgs->J, &J_local));
961       PetscCall(MatSolve(J_local, lbfgs->rwork2_local, lbfgs->rwork3_local));
963     PetscCall(VecRestoreLocalVector(lbfgs->rwork3, lbfgs->rwork3_local));
964     PetscCall(VecRestoreLocalVectorRead(lbfgs->rwork2, lbfgs->rwork2_local));
965     PetscCall(VecScale(lbfgs->rwork3, -1.0));
967     PetscCall(MatMult(lbfgs->YtS_triu_strict, lbfgs->rwork3, lbfgs->rwork2));
968     PetscCall(VecPointwiseMult(lbfgs->rwork2, lbfgs->rwork2, lbfgs->inv_diag_vec));
969     PetscCall(VecAXPY(lbfgs->rwork1, 1.0, lbfgs->rwork2));
971     if (lbfgs->strategy == MAT_LMVM_DENSE_REORDER) {
972       PetscCall(VecHistoryOrderToRecycleOrder(B, lbfgs->rwork1, lbfgs->num_updates, lbfgs->cyclic_work_vec));
973       PetscCall(VecHistoryOrderToRecycleOrder(B, lbfgs->rwork3, lbfgs->num_updates, lbfgs->cyclic_work_vec));
976     PetscCall(MatMultAddColumnRange(Yfull, lbfgs->rwork1, Z, Z, 0, h));
977     lbfgs->Y_count++;
978     PetscCall(MatMultAddColumnRange(lbfgs->BS, lbfgs->rwork3, Z, Z, 0, h));
979     lbfgs->S_count++;
992   Mat_DQN  *lbfgs;
1015   PetscCall(PetscNew(&lbfgs));
1016   lmvm->ctx              = (void *)lbfgs;
1017   lbfgs->allocated       = PETSC_FALSE;
1018   lbfgs->use_recursive   = PETSC_TRUE;
1019   lbfgs->needPQ          = PETSC_TRUE;
1020   lbfgs->watchdog        = 0;
1021   lbfgs->max_seq_rejects = lmvm->m / 2;
1022   lbfgs->strategy        = MAT_LMVM_DENSE_INPLACE;
1024   PetscCall(SymBroydenRescaleCreate(&lbfgs->rescale));